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considered include: hydrodynamic systems; chemical, systems ana 
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structural mechanics; njDndestructive testing and tomographic 
reconstruction; mathematical models in the biological sciences; ,and 
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O^RVEW 

The 1980s are a time of profound challenge to the t6chnolo|ical strength ^ * 
of the United States in the economic as well as the military spheres, and , 
,our countr> 's performance in research and development in its engineering , 
laboratorie*^, be an important and perhaps determining aspect of our 
success in meeting this challenge. Advanced engineering develop^nent is 
now based mainly on scientific computation, which in turn relies on math^ 
ematical modeling and laboratory experiments,^ Together th^y represent 
one of the areas in which the strength of^ nations is being tested today. 

Mathematics is essential in the development of theoretical and com- 
putational models for solving the highly complex.pioblems of engineering 
and basic science, which encompass a range of scfentific difficulties^ On one 
side are questions of computer architecture and the science of aJgorithms. 
. On the other Side is the modeling of chemical and physical processes 
* by means of ihatheinatjcal equations. The issues are tied together by 
mathematical theory, ^hich seeks a fulTunderstaiKlitig of the nonlinear 
phenomena eontained in the equations and implements ttfis understanding 
throogh computational methods. This span of sci^tific activitiesjorms 
the subject that is known as applied mathennttics. 

The pilncipal conclusion of this committee is that computational 
modeling, which is a hi^-ieverage element jof our nation's scientific and . 
technological effort^ requires increased emphasis and support. The conclu- 
*'sion is documented by an examination of typicakapplication areas, which 
reveals the pervasive difficuUies that accompany computation of realis- 
tic problems and leads one to consider both what computers can do and 
what they cannot currently do but might eventually be capable of doing, 
As illustration, we examine several <3eep theoretical problems, iiicluding 
turbulence and combustion. At the frontiers of attack on these prob- 
lems we discover the limitations imposed by our current understanding of^ 
model formulation and computational capability. We examine modeling 
^problems and dgorithms a^iapted both specific features of the desired 
solution and to th^ computer architecture. We also examine compuler 
architecture alid component de^gn and manufacture as a mathematical 
modeling problem. 

The Committee recommends increased support for 

1 Hesearcb'in compjutationat modeling and applied mathematics, 
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■ ' 2, Cpmputmg facilities dedicated to this area^ . 

3, -Education and manpower development in computational and ap- 
plied mathematics. ' . ^ \ 

These recomimended increases include fli>ancia) support from gojern- 
ment and industr>^ institutionat support from universities. * ^ 
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l;NTRODUCtION 



The extenyue use of compu,ta'S in advanced development work began dur- 
ing World War and today computing is a vitak component of science^ 
engineering, and modern technology. Most advanced technological devel- 
opment, if rom aircraft design to automobiles to petroleum to satellites now 
follows this pattern of reliance' on thje computer. Moreover, the needs of 
national defense have posed scientific and engineering design problems as, 
difficult as an>*ever encountered. Numerical computation and applied 
" mathematics have played ^ essential role in dealing with suci) problems. 
In fact, numerical fluid dynamites was bom during the 1940? for the pur- 
pose of assbting iif the design of nuclear iveapons^ combat aircraft, and 
conventional ordnance aiiti is now applied widely by industry^ 

Most problems of engineering or scientific interest are too complex 
to be modeled and compute;! exactly* Instead one considers a series 
of approximate models' and computations, each of which illuminates a* 
different aspect or idealized portion of the overall problem. When used by 
a skilled engineer or scientist, these mathematic^lmodels greatly enhance 
tlxe judgment that goes into design 'decisions ao^educe the amount of 
expensive laboratory^ and field testing required. These advantages account 
for the widespread use bf these models, ^ * 

More specifically, mathematical models are used in engineering design 
problems in the tollowing modes: ^ ' . / 

h To provide the first tQst of a new de^^ijdea. Beyond common 
sense and simple hand calculatiops, the computer model is usually the^ 
chef^ltest and fastest test of an idea* Thi^ test is applied before deciding 
whether to conduct a series of experiments or to builds prototype. 

2. To reduce the time and coat associated with laboratory and 
field tests. Usuajly engineering problems contain several critical des^ ^ 
parameters that will have involved a certain amount C(f trial efid error in 
the search for the optimum choice. The B5mputer is used qualitatively 
(Will an increase in parameter X improve or degrade performance - or ^ 
performance parameter y?) and quantitatively (Whicl^values of Klesign 
parameters JCi , . . JCn will op^imiz^ performance?). 

3. To assist*^in laboratory or field tests that determine model 
parameters equations. Often the model parameters are measured 



only indirectly. Thus, a coiApul^r' model of the laSoratory apparatus 
may be needed to extract the desired information from the observed data. 
Usually these models are simpler than the complex engineering models^ 
and their defining equations c^n be solved with greater precision and fewer 
approximations. 

4, To replace laboratory or fi^ tests/ Sometime^ tests are impossible 
or Impracticable. For example, measurement of chemical reaction rates 
' ^t extreme conditions of pr^sure and temperature is Very dif{icuU» and 
accumulated experience through trial an9 error is not adequate for soKJng 
the probrem of landing the Space Shuttle, for reasons .6f human safety as 
/ , well as cost, ^ * f * . 

5 To improve the education and Judgment of engineers and 
tists using the models. The mathematical models^d {Computer solutms 
^ provide a vast increase^ in the quality and caliber of the data. Thus while 
the laboratory measurement may produce some overall quantit> (e.g., total, 
flow in and out), the computer model might yield detailed velocities and 
concentrations at each^jooint of the flow field* Because the equ^itions are 
nonlinear, it is^ifficult to foresee all the relevant phenomena, much less to 
uhderstand their relative importance. Thus the mc^del becomes for u^r$ 
an Experimental tool that allows them to understand a problem at a level 
of^detait that cannot be achieved by other meatvs. 



In summary, mathematical and computer models arfe^ised because 
they.are faster, cheaper, and more -effective. ,^ * 

However, tnodels have limitations. There are limitations in the 
validity of the equations used, in the adequacy of the solution algorithm, 
and in the si^e and speed of the computer, AIso^ the cost, accessibility, aUd 
reliability of computer softflyare and, sometime]^ the cost of the computa- 
tion itself can be liiriiting factors. These limitations in some sense define 
'the frontiers of science, but more specifically they define the frontiers of 
applied mathematical science, 

Problerils of realistic interest typically involve the study of diverse 
physical phenomena ori many scales of length in fully three-dimensional 
settings. Though essential, experimental science in these contexts i^ ex- 
pensive and difhcutt. The design of modern strategic weapons $ystems 
epitomises these characte'ristics. For example, in the de^gn of a TRIDENT 
submarine, the architecture of the vessel, tl^e design -of -the *missite, and 
the design of the nuclear j^rathead all need t<i t)e iripdfeled and integrated, 
which Quires many tKoujsands of hours of c'ts^puler time and stretchy 
available copiputing power and modeling techniquesHo their limits. 

Generally,. desim and jevaluation of new kinds ofl^defeuse^systems^ 
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$uch as remote!}^ controller} or robotic ve^cles/ ivil] require analysis of 
entirel> no\el ^^tems, using paraJlel computer architectures among other 
things. Tlmn^ds of '|ndusti;y for tethn^logical advan<;es are similar and 
are dommaied b> suclv basic concerns as poUu^ion, depletion of resources, 
eogrgy conservation, and efficient use of manpower. It seems safe to 
say that defense and industrial needs will continue to lead numer|:al 
computation 'and applied mathematics into new and chaltetiging regimes 
Depending on the complexity oftheptoblem aird the magnitude of the 
effort exptnded> models range fropi excellent to merely suggestive in their 
qualit> and usefulness. Jn all cases, improvement of completer, modeling is 
joxi^ of the most promising avenues to improved techtiol^ca! performance 
^ by our nation. . 

■|L As one of the aims oLthis report, the Committee wantBL^ show and 
^ emphasize that in the computational approaches to most of today's press- 
< Jng and challenging scientific and technological problems the mathematical 
aspects cc^ot and should not be considered in isolation. Th^re is a unity 
' aqiong the various steps of the overall modeling proces^rom the formula- 
« Uonof the phy^al problem io the construction of appropriate' mathemati- 
cal models, tlf44sMg^ of suitable numerical methods,.tHeir computational 
]mplementection> and, last but not least, the validation and interpretation 
of the computed results. In particular, the Committee wants to illustrate ' 
that the step| are more often than not deeply interconnected and that 
the computational process may iniil^d be pa^^Df'^ie model construction. 
At the sam^time, there are problem areas, such as turbulence, where cur- 
rent theoretical research- may promise a deeper insight into an important 
physical phenomenon. A . - * , , 

In litie with these^ms, the report uses a '^matrix approach^* that 
views the ^ame problem, from three different standpoints. In Chapter 2, the 
traditional approach is taken/of discussing a number of typical problems 
' . leading^to computational modding from the viewpoint of the acientiflc 
or engineering discipline in which they arise. Then in Chapter 3 certain 
^ of these^problems are touched on again, this time from the viewpoint of ^ 
the^ computational and mathematical di0ieulties that arise in connection 
with them. For example, these difficulties meiy involve latge number^ of 
degrees of freedom; different scal^ of time and length, or "singularities of ' 
various -types. In Chapter 4 the viewpoint becomes that of the numerical 
^ ' algorithms involved in th^ computations, .such as various discretizSttm^||^' 
inetho^ j<$htinuat]on approaches, and splitting tecbiuques.. 

CSjiecessity, many topics have been left out. The list of applications, ■ ^ 
fpr example, is by no means complete, ^nd, in fact, entire ar^as such 
as.reactor safety and reactor physics are not mentioned at all. Neither 
did the.Committee attempt address all com^tutational and matheinati- 
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c^l dilBculttes ^br all variaUons of numericat algorithms. The choice gf^- 
|opics Avas motivated in part^l^y erxergy-r^ted considerations, the exper- 
tise aiid interests of the Committee and its^visers^ and. the report's broad/ 
purpose, Avhos^ achievement would be hindejred by my attempt to be en- 
cyclopedic. For the same reason, the repott is ndt;;itrfflEIftkc|^ be a t^h- 
nical summaiy^ and this is als^ reflected in tl^ fact that no attempt was 
madejkt reference the relevant literature. The repoH is mainly addressed 
to sci^iifically literate Yeaders who know how to consult the literature 
whenr-necesary. ■» ' ' " 

' "^The Committee obtained advice and technical support from many 
colle|Lgues ach^ss the countiy and from abroad. A li^ of names of^^aJl 
'ttio^who helped in this Y^ork Is given In Appendix A, and the Committee 
is esttremely grateful for all th^ often extensiv^ documentations^ special 
write-ups, and other comments thaC ^.received* The report is the r^ult of 
a study he^n by the C6mmi:^t!^!in 1981 on cbmputational mathematical 
modeling ahdfnath^matics applied t^the physical sciences with particular 
reference tb tji^eads of tfie Departrfient Energy (DOE). The preparar 
tion of th^ l^e^rt wa^T^port^ii by the Applied Mathematical Sciences 
Research Program^ of the Office of Basic Energy Sciences of DOE, and the 
Committee aJso expresses its thanks and appreciation for Qils support. 

As stressed ii^ th^Overview, the Committee found that improvement 
of the m^thepia^tcal and' Computer modeling of scientific problems ia an 
• important priority for our nation. The cl^lenge is broad, and there are 
i;o simple remedies for current shortcomings. Accordingly, the Committee 
recommends the following: 

1 ^ftthased reseaf^^uppoTt for computational modebng and applied 
mathematics . . 

The technological challenges of the coming decades will impose new 
tests of our abilities in computational and applied mathematics, and meet^ 
ing the tests will require increased research effort. M illustrate in this 
report, the challenges are typically multidiscipUnary in nature with applied 
mathematics and modeling often in a central position. ; * v 

H^nc^, to support research in this area, multidisciplinary teams of 
adequate size to make progress onthfi^fi'^eGimpl^^ problems should be 
encouraged; and olrganizational means shout^e devised to facilitate their 
establishment^ continuity, and success.. 




5. Increased support for computing facilities dedicated to compiHa- 
tional modeling and applied mathematics ' 

'Ready access to- modem ccmpuiing systems is essential. Laitcjjf 
equipment is the critical jW:tor.most strongly limiting academic researct^ 
in computation^ mathem^cs. There is need both for conveniently usable 
local equipment and for access to large-scale computers. Local facilltt^ 
are necessary for entire problems of modest size and for such taso as 
code development, interactive debuggings test runs^ andtgraphical analysis. 
Large-scale computing is essential because of the size of many ^ical 
problema as documented iiuthis report. In this connection^ the Committee 
strongl} endorses the findings ^d recommendations in the recent Report, 
of the Panel on Large-Scak Computing in Science and Engineering, Peter 
1). LaXf Chairman^ National Science Foundation, December 26, 1982, 

> 

5, Increased support for education and manpower developmerU 
computattond afid applied fnathematics ' ^ 

^ Todays there are unmet manpower needs in computational and ap- 
plied mathematics, as discussed^ for example, in Science and Engineering 
Education for the 80s and Beyoy^ a National Science Foundation report^ 
October 1980. These needs are foun^ in industry^ government laboratories^ 
and academic institutions. The cntical challenges in this area call for a 
focus oh quality. Specifica&y^ the complex interdiacipfinary nature of the 
problems poses special educational challenges for students and young re* 
searchers, and graduate and postdoctoral fdlowshlp support for participa- 
tion in multidisciplinary teams of the ^e.discussed above woiild also be 
helgtful 

A valuable aspect of such multidisciplinary education is the/mterac- 
tion it creates between applied mathematicians and other applied scientists 
in universities, government, and industry. 



2. APPLICATIONS 
2.1 HYDRODYNAMIC SYSTEMS ■ - 

Hydrodynamic proceases touch nearly every aspect of our lives^. In a report 
even many times larger than the present one^ we coujdrtot possibly discuss 
all hydrodyBamics api5licatioDs, Therefore, we have chosen a few^ with 
the hope that they will serve to illustrate the importance of hydrodynamic 
models and to point out some of the projiilems that occur m developing 
Ahe models. , 

A large fraction of the lurrent research in computational modeling 
iadone with hydrodynamics applications in mind. The great variety of 
E^ponses that we received to our requests for material for this report 
^test to the diversity of tfies^ appli(^tions» Among them were applications 
haying to do with aircraft and wing design^ both at subsonic and supersonic 
speeds; global weather prediction and local phenomena such as tornadoes^ 
water wav^s and ship bull design; piping networks, such as in nuclear 
reactor or power plant design; geologic phenomena, such as glacier Sow 
or convection in the Earih^^s mantle; biological flows^ such as the ficfw of 
blood in the heart; and chemically reacting flow, suck as combustion. 

The general system of equations governing hydrodynamics are called 
the Naviei^$tokes equations. Theyare.a stateinent of mass and momen- 
tum conservation^ the momentum equation being a fonnuiation of 
Neon's second law, F - ma. The Navier-Stokes equations were first 
developed in Prance by Navier in the early iSOOs. They represented an 
improvement over the Euier equations that were first derived in 1755, in 
that the Navier-Stokes equatims included viscpus effects that were absent 
in the Euler equations* However, it was not until 1904, when P^^dtl de- 
veloped the botmdary-layer approximations, that predictions of practical 
viscous flows coi^ld be made. Practical solution of the f^ Navier-Stokes 
equations had tp amit the development of modem bigh'speed compliters 
beginntiig in the 1940s. ^ ■ / . 

Why are equations that have been known for over 100 years so hard to 
solve? The answer lies largely in their inherent nonlinear characteristics. 
]ii!unediately upon looking at the equations one sees that the convective 
transport terms (the acceleration in F = ma) involve velocity times its 

V 
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gradient. This nonlinearity is always present/ and it is respojisible (or 
the existence of complex phenomena sucb as shoek waves and turbulence. 
In principle the Ni^vier-Stokes equations aloni provide a description of 
turbulence, however, one would have to resolye.such small length-scales in 
their solution that this appro^h is not pf practiod importance. Therefore, 
a gr^at many approaches to approximating Wbulehce ejects are being 
pursti^. Typically, these models introdiice further nonlipearities into the 
system. 

Other problems that arise in the solution of hydrodynamic problems 
are related to the disparate time and length scales that must be resolved. 
In particular, convective transport is characterized hy the fluid velocity, 
whereas pressure waves travel at the sound speed-typically orders of mag- 
nitude faster tha^uid speeds. At the same time the effects of diiiusive 
processes {e.g., stmring stresses] are felt instantaneously throughout tb? 
flow. In some cases the fluid can react <;hemically.. In those cases^ the reac- 
tion rates display a strong nonlinear dependence on the fluid temperature. 
This Introduces yet more characteristic scales into the models. 

Depending on characteristic parameters such as the Reynolds number, 
the solutions to the Navier-Stokes equations either can be smooth and 
steady or they can exhibit regular oscillations. Or they can be completely 
chaotic. It is clear that depending on the regime the appropriate solution 
procedures could be quite different. 

Finally a word about totodary conditions. In some proh^s the flow 
is enclosed, and hence boundary conditions are applied at the boundaries 
of^he enclosure (e.g., an automobile engine cylinder), which, can often 
have a complex shape. In other cases^ such as an aircraft wing, the flow 
is effectively unbounded, and the boundary condition should be applied at 
"infinity" (see Section 3.4). Some modeling problems arise in a^proximat- 
, ing infinity by some finite boundary. In still other problems, such as 
the flow of blood in a heart, the boundary is both complex in shape and 
deforms depending on the forces e:!certed on it by the fluid. Modeling such 
a problem is clearly a challenga 



24,1 Wings, Wind T\innels, and Computers 



The economics pf the energy shortage implies that planes will fly at speeds 
close to hut less than the speed of sound, At such si>eeds there , 
**supersonic bubBF*^" over tlje wing, where the local velocity of the^ 
relative to the wing is greater than the ,3peed of sound. In this case the 
presence of shock waves is typical and undesirable. They are undesirable 



■ (aV •. • (b) ^ ■ ; (c) V' . ^ 

FIGURE 24 Distribution of velocity (normalized) on surface of a Korn supercritical airfoil at Mach numbers 0.750, 
0.752, and 0,748. The long bar at?i^ is the transition value. {PYom C, S- Morawetz, The mathematical approach to 
the sonic barrier^ BvlL Am. Maih. Soc. 6, 1275^143 (1882)-] 



because the drag can be compi^ted as being proportional to the third power 
of the $hock strength. The gosl of efiicieat wing design is to produce 
wing shapes with no shocks or only weak shocks in thi$ transonic region. 
A general mathematical theory shows that shockle^ss wing foils exist for 
given transonic cruising' speeds.. However, the problem of finding such wiag 
shapes is both overdetermined and extremely sensitive to small changes in 
the data, i,e,, "ill-posed" (see Section 3*6)* The solutitih to such ill-posed 
problems is still valid from an engineering point of view because operation 
at neighboring off-design conditions will produce only weak shocks and 
small drag, see Fig. 2.L . ■ ' - 

Using computer algorithms created by applied mathematicians, it is 
now possible to solve both the inverse problem (design) and the forward 
problem (of determining the flow field for a given wing shape and velocity) 
lyith sufficient accuracy that the use of costly .wind-tunnel experiments can 
be greatly reduced. This accomplishment is a striking success of recent 
appiied mathemi^tics. / 

the theoretical areas that have contributed to this study include the 
theory of nonlinear elliptic equations, complex function theory, andmixed 
problems, for which a prototype is the TWcomi equation, , 

which is elliptic for > 0 and hyperbolic for j^<0. The ellipticgsregion 
' corresponds to the subsonic region, and the hyperbolic region to the su- 
personic "bubble/' In a numerical method for the design problem^ an 
analytic continuation makes the elliptic region: hyperbolic. The resulting 
equation is solved .by . the method of characteristics, and then the analytic 
contiiuiation back to real values of space is performed numerically- 



Chatis, l>irbulence, and Droplet ' \ ^ ' - 

- - ^ ■ ' -s^\. ^ ^ " : ■ ^ ' ' 

, Turbulence produces a boundary layer along the trailing edge of an aircraft 
i^ng. The boundary layer degrades the wing^perfjcjmance and thus is,an 
importapt part gf t,)ie.d^i^ protilem. The flow behind the trailing edge of 
a wing contains a vortex:""§hee^j and th? roll-up of this yortex sheet produces 
tu^bulence.that confttitufea a.safety hazard for small aircraft flying in the 
,wake of a jumbi) jet». Tfie axi^ptlhe yortex rojU-up Is perpendicular to the 
wing, ^d so the rollrup is intrinsically thr^e^'dimension^. tsi the simpler 
case of twoidnnensional .flow the vortex ;sheet is a line or curve in. the 
plan.e^ emanating^from the trailing edge of the wing foil. In r^ns in 

- . ^ , . . Jl " . . 
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which th^ line is streUhing, it is geometrically atable. InTegions in which 
it is contracting, it is unstable. Instead of contracting^ it fon^s a spiral- 
vortex structure and hence is stable. 

The two-phase floW of water and steam in a cooling pipe, or of oil 
and gas in an oil-^reservoir production well, is also a problem in wbich 
the geometrical instabilities of large-scale fluid motion are important, 
Here,an internal movab|p boundaiy separates regions of different material 
■ properties. In some cases (a heavy Suid, e.g,, water, falling into a light 
fluid, e.g., air^ the material interfa^eJs unstable against formation of 
fingers. There is continiiation of the nonlinear Instabilities leading to 
pinchqfT of droplets and a chaotic regime (mist) that can be analyzed on 
various length scales, as discussed below in the case of turbulence, ^ 

During the combustion stroke of an automobile engine, the flame 
is quenched^^when it reaches the cold cylinder^-walls, and incompletely 
burned/uel present in the combustion chamber at this tipe contributes to 
rjrdlution ai>d to a loss of fuel ^ciency. The rate of advance of the flame 
' from) the spark plug to the cylinder walls is governed by the laminar flame 
speed and the rate of turtjulent mixing. Of these two effects, the second 
, is more important. The turbulent mixing is produced by vortices that 
det^h from tlie turbulent Jboundaty layer at the wall during the intake 
and compression strokes (see Fig, ^ 2.2), Thus, an accurate 'modeling of. 
this problem requires an ability to treat a number of fluid Angularities: 
flame fronts, vortices, turbulence, boundaty layers, and>J>ouiidaiy<*l83rer 
separation, * ^-^j, 

The examples above dhow that sing\4arities in fluid flow' may be 
geometrically unstable. When this Instability occurs in a regime governed 
by the scale-invariant Euler equations^ the phenomenon is repeated on 
all length scales and leads to chaotic ^ohitions, IVirbulence^^ vortex roll* 
up, convection finge^ringj and droplet formation are examples 6f thi^ 
phenomenon, which we now discuss from a'^eneral point of view. 

The Euler equations of fluid dynamicjf allow intrinsic singularities, 
namely vortices, boundaty and shear discontinuity layersj i^ontact or 
material interface discontinuities, and shock ^aves (Fig, 2,3 shows a com- 
putation of the stretching of a vortex tube in a periodic inviscid flow). 
Depending on the problem, special discontinuities such as flame fronts or 
chemical reaction fronts (for fluid dynamics with chemistiy) may occur. 
Within the singularity, the Euler e<^ations fail to be a correct description 
of nature, and corrections (either parabolic effects or perhaps a more com* 
plicated Euler equation with more state variables) may.be required,. As an 
example, consider a shear l^r (i,e,, a jump discontinuity in the tangential 
velocity component). Taking the curl of the Navier^Stokes equations, we 
. obtain for c ' 
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w= Vxv = vortidty 



the equation 

" *'^ + (V'V)(j-(w'V)v = i/Aa;* ' ' 

where 1/ is the kinematic viscosity. To understand the significance of this 
equation^ we specialize to two dimensions. Then w is a vector in the z 
directi(yi, {uj - V)v, and*' 

is the total, or Lagrangisn time derivative, so that the Navier^Stokee 
equatioii^ says that yotticity moves by passive transport plus diffusion. 
The extta term, (w - V)v above, induces vortex production as a result of 
the stretching of vortex lines in three dimensionar'and is important for 
considerations of geometrical stability as discussed below. In. summaiy, 
the diffusion, term i/A^ of the Navier-Stokes equation is responsible for 
the vorticity leaving a l>ouxidaiy or internal eh^ leyer and diffhsing mto 
the rest of the flow- Without viscosity there is no mechanism for voi^city 
to enter (an initi^ly irrotational) flow region. The Prandtl Iwundary-leyer ^ 




FIGURE 2.2 Stretching of a flame by avgrtical structure* Such stretcliing 
is important for the efScient operation of.engin^s; it enlia|iees,^|^giiing by 
increasing the ar^ of theJame. (F^om A, J. Chorin^ Flame ad^ni^tion and 
pppagation algorithms, J, Cqmput Pkys, 35, I'll (IJBO).) 
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FIGURE 2,3, Successive stages in the stretching of, a vortex tube in an 
inviscid periodic flow. Vortex stretching is the mechanism by which energy 
in a turbulent flow is transferred to eyer smaller scales Where it is eventually 
dissipated, a, step 10, tinle = Q.65; % step 20, time ^10.88; c,,step j30, 
time \ L04; jd, step 40, time — I-2L (FVom A, JJ'Chorin, The evolution 
Qf a tufbulent vortex, Commun, Math, P/tvs,'83,, 526-527 (r982).] ' 



equations are a special version of the Navier-Siokes equations (scaled in the 
4iortnai direction, s^ that, diffusion Ojccurs only normal to the boundary). 

Often the fluid singularities are ^eometriciflly unstable. They may 
bifurcate in a predictable fashion, developing ^^rolls*" ^ouett^.flow) or 
"cf/ils^(^nard flow), or they may become irregular and higbl^ convoluted 
with a tendency toward chaos^ known as turbulence. There is no scientific 
reason to question the^idity of the Navier-'Stokes equation as a micro^ 
scOpic description of physics even into the turbulent regime. However, its 
usefulness^ a description of large-scale fluid motions in the turbulent 
region can W questioned, and some other description of fluid Sow, such as 
a random ensemble pf interacting vortices, c^Sul^ be more effective, 

there are^ree enjergy orlen^h sdftlesTB which quite distinct charac- 
teristic phenomena dominate, T?he smallest length scale is that on which/ 
energy dissipation dominates. On this length scale, the Navier-Stokes 
equations are ih^ correct equations. The viscosity is large^ causing velocity 
fluctuations tp be fapidly smoothed and solutions to be (locally) ^laminar,^* 
Since '^oothness of solutiops is a local. property, we may conjecture that 
alfsolutipns of the Nay^er-Siokes ,equatioi^ should be smooth for^l time. 
This conjecture is the major unsolved probl^ of the energy dissipation 
r^ge; Jt U krbwn that solutions^ with smooth data will remain smooth 
for a short time, and si^lutipns with smooth ai)d small data will remain 
smooth for all tiihe. Both statements exclude turbulent regimes. 

For weak solutions, Leray showed that the^set of times t Sot which 
v{xjt) fails to ^^mpoth i^ a ^t of theasure 2ero! Recently-considerable 
progress has boen^made in r^trictihg the possible singularities of the 
jNavier-Stokes equations, - v ^ - , 

F^r longer length scales, viscous effects do not nlay a^ajor role 
and the fluid' flow qan described by Euler equationVN^ovrever, this. 
3implificatio|i gives rise to problems. The problems are notrmerel^r techni- 
cal but reflect the intrinsically complex phenomenology of nuid dynamics^ 
The Euler equations are scale-invariant. Thus if v v(x, t) is a solution, so 
is V » v(a:r,a£). The inertial range is the's^t of length scales dominated by 
scale-invafiant, uniWsal physics. Whatever phenomena can occur (e.g,, 
vortices) will be repeated on length scales in the^'inertial range. The 
inertial ran^is limited at the smalkr^nd^by viscous dissipation. A) the 
larger end, Wis limi^^ by'the special bdiindary joad, initial coilditions 
imposed on the flp^w, which result In special' (nonun{versal),flow behavior. 
' * Thejnertial region is dominated by ^ale-mvariant behavior, ^here 
is a flow of energy from large-scale motions to smalleiv onea (the **eoer^ 
cascade"/. This cascade seems plausible on physical grounds as a of 
third law of thermodynamics hut does not have a . rigorous mathematical 
status. It can be explained as a consequence of the^eometrical instability* 
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of vortex lines and jshear layers. As these go" unstable, they generate 
(smUl^r) new vortices" and vorticitjt- 

The energy cascade leads to a dtimensional analysis o^characteristic 
exponents and^to the Kotnu)gorofl[5/3 power law 



for the energy distribution as ^function of mve number k. discussion of 
the experimental data in connection with the Kolmogoroff theory has the 
vortices, which occur on all length scales in tlie inertia! range, as* filling 
space. Actually, it may ]be better to assume the co'htrary:," vortices of 
a given size fill only a small part of space. Then the sinaller vertices, 
which are driven by the largenone5, will occur only within the region of* 
these larger vortices, and in fact only within Ismail part of this region. 
This is the notion of^tennittency. It leads to the idea of a singular set 
for solutions of the EjjJ^r equations, which is a Cantpr set of fractional 
dimension less thaifX^* , ^ ^ * , 

' Intermittency leads to modifications in the Kolmd^rofF exponent 
and to a renormalization group type picture of turbulence. Numerical 
calculations to determine intermittency and energy cascade .exponents ' 
have been performed. The calculations start by tracking vortex lines in 
a fluid flo^^nd proceedJJux>ugh a sequence or renormalization group 
length scale transformations Jo focus pn the singular Cantor set within 
tEe solution. 

In most proWemSt the inertial region contains lei^hj too small tor be - 
'used directly in a fluid calculatioir. Its importance lies in its role of fixing 
parameters feuch an effective or eddy viscosity, which are then used to 
determine the large-scale motion of the fluid. The inertial region is not 
particularly well understood from either the theoretical or numerical point 
of view. ' 



The large-scale fluid motions are produced directly by the iniUat and 
6oundary^j:on(litioAs imposed On the flow. These motions are strongly 
problem dependent. Numerical calculations and esqjeriment are important 
tools in their study as is the detailed analysis of simplified and idealised 
flow configurations. An important theoretical question is the evolution of 
initially unstable flow configurations. This qu^ion arises in connection 
^Yith the pnset of turbulence and in connection with the energy cascade, 
where large-scale vortices excite and drive small-scale ones. ^ 

4n some problems (supercritical turbulence), the instability in an ini- 
^ * tially laminar flow is nonturbujent but arises f^n^ %e bifurcation of a 
fixed point FVirther bifurcations lead to^bigher^dimensional torit and a 
general theory explainslthat g^nerically the Sow on the (sufficiently high- 
dimensional) torus has, a strange attractor as its limit set and that this, 

^ ■ m 
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strange Sttractor is chaotic in nature. This picture has been analyzed in 
the context of the Lorentz flow, which is the tnincatjpn of the Navier- 
Stokfes equations to include only a small number of m^des,^ The strange 
attractors found there have a Cantorlike structure. An example of super- ^ 
critical turbulence is Couette flow. ' 

Subcriticai turbulence occurs when the finite (noninfinitestmal) per* ' 
turbation b less stable than the infinitesimal one. Then turbulence occurs 
below t^ie critical Reynolds numbers at which the linear theory' shows 
-instability and may go directly to turbulent behavior without a discrete 
sequence of nonturbulent bifurcations starting with laminar flow. 

2.2 CHEMICAL SYSTEMS AM3 -COMBUSTION i 

FVom thet invention and manufacture of an enormous range of syn/ 
thetic materials (e.g., pfastics) to the refining ^d burning of fossil fuelsA 
chemistry and ch€fnicaJ processes affect nearly every phase ot our liyes. 
Naturally it is important to understand ^d control, as fully as possible) 
many of these complex chemical processes. For example, we seek to find 
new and ^tter materials, to reduce costs^ and to generate energy more^ 
efficiently and with less pollution. Appl^ mathematics and computa^ 
tional modeling d)ntinue to play a valuabft role in meeting these goals. . 

One of the oldest chemicarprocesses harnessed by man is combustion. 
The* successful modeling of combustion provides an extraordinarily rich 
.source of challenges for the computational mathematician. Frequently, 
combustion models have to incorporate all the difficulties of compUcated 
fluid me<^Shics coupled with complex chemical kinetics. The challen^ 
include developing algorithm^49 ^^^^^ accuracy and to reduce computer 
time and stoi:^e. The modeler^tlso seeks appropriate simplffications that 
take advantage of any special attributes of the governing physics in order 
to gain more efficient computation. Since combustion contains a wide 
range'of chemical processes, we use it here a| an example to illustrate 
points of mathematical interest in general'chemicaJ systems. 

Even within ^the topic of combustion there is an enormous diversity 
of applications. The first topic that probably comes to mind is the mode^ 
ing of internal combqstioii engines, and this is an impQrtaot. application. 
Modeling is a^art .of ongoing research to design new types of engines (e.g., 
direct injected stratified charge)) to improve fuel economy, to utilize alter 
nate fuels (e.g., alcohols), ^aud to reduce pollutant formation. Sifnilarly, 
research for other combustora, such as gas turbines or power-plant boilers 

. "it . . ^ 
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benef|,t Trom computational models.. Still there are maii> othenrimportant 
combustion problems Aide from power generatioR. For exaiflple^ the field 
of fire research is devoted to problems'such a^ bonfires spread in^buitd- 
ings and the behavior o^atious fire-retard$jit materials. An imfportant 
current topic in reactor sai^y is the characteriz^ion of hydrogen-air fires 
such as occurred in the Thre^ Mile Island accident. *\nother pubKc safety 
question deals ^ith the problems of fire atud explosion associated with 
a liqgid-nSTUral-gas tanker accident/^roblems of burning Poal an^j coal 
gasification are also Eopics of great current interest. ^ ^ v 

Perhaps simplest chemical process tronS a physical viewpoint is 
chemical equilibriilhi. At equifibriuhi all chemical reactions are a&sumejd 
to have gone to c<lnipletiont and the specie? concentrations ^are such that 
the mixture is in a minimum free-energy state. The mathematical com- 
putation of the chemical equilibrium state is posed as a coi^straitied piinl* 
muatioa problem. In combustion the eqijilibijum composition corresponds 
to the products of combustion long after the combustion is complete. 
, Phjsicatt]^, the next step in complication comes ^^ith the inclusion of finite- 
rate chemical reactions. Here the mathematical j)roblem is one of solving 
systems of stiff ordinary differential equations (see Section 3.2 for a discus- 
sion of stiffness). Models^ of shock tube^ or ffow reactor^^ which ar^ used 
frequently to probe fundamental questions in chemical reaction kinetics, 
fall into this class of problems. The phjrsical situation is complicated fur 
ther b> tiie iticlusion of fluid motion and heat and mass transport. In this 
case the mathematical problem isV)ne of solving systems of p&rabolic*or 
elliptic partial differential equation's. 

Consider the intemar combustion engine as an example. What are 
the things that we miglit hope to learn from modeling? Ultimately, we 
hope to influence geometrical considerations such as^ combustion chamber 
shape and component placement (e.g^ valves^ sp^tk>BtugSf fuel injectors). 
vie also hope to pi^vide a fundament^ understandinki on the molecular 
level, about how fuels are oxidized and how pollutants are formed. With 
such understaA^ing we caiv suggest ways to alter the ccjmbu^tion proce^ 
to advantage. In the past only power and size were important considera* 
tions, and engine optimization could proceed mostly experimentally. Now^ 
however^ there are too many parameters to optimize simulj^eously, and 
computational modeling is increasingly important. 

Because of limitations in available computational resources, two tacks 
^ taken in engine modeling. One is to consider mostly ^bydrodynamic 
effects. Here the modeling of boundary shapes and component placement 
is of ^primary importance (e.g„ How should the piston face be aji^ped, 
and where should the SQM*k plug be placed?]. Coftiplex domalna and two 
and three-dimensional enects are important. The models must incorporate 
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-FIGURE 14 Velocity vectors computed in a direcMnjection, stratified- 
charge engine at a position near top dead center. The combustipn takes 
place in a swirling environment in a cup-like r^an machined into the 
pisU/h^ Several vortices%^ ^eea to develop in the cup* [From T- D. Butler, 
L- D. Cloutmaft, J. K* Dukowicz, and J. D^ Kamshaw, ^Juitidimeasional 
numerical simulation of reactive ii^intemal combustion engines, Proj. 
Entm Combust Sci\^%3^ ti?S?_ffl/ i^^^ ' * 



turbulent hydrod^amic ef]fecis\,^d^ sometimes phase-change processes, 
such as fuel spfiay inj^tions. The chemistry is ueu&lly simplified in these 
models becauseit is not feasible to consider both complex chemical kinet^ics 
and hydrodynamics on current computers. Figure 2A shows the velocity 
vectors that result from a two-dimensional simulation of a direct injected, 
stratified charge engine. This i&^ new engine concept that is being studied 
in a U*S- Department of Energy-sponsored cooperative pro-am including 
General M9tors Research Laboratories, Princeton University, and three" 
National Laboratories. 
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In addition to the h3^rodynaniic aspects of engine combustion, there 
are impoTtant unanswered questions about the chemistry. Therefore, the 
second tack is to consider simplified hydrodyiianuc situations^ such as 
lanUnar flames, and treat the chemical kinetics in great detail. 'These 
moctels address issues such as ignition p)ienomena and pollutant formation^; . 
Figure 2.5 shows some species profiles computed in an atmospheric pres- 
sure acetylene-oxygpn prembced flame (acetylepe is an important reactant 
\n soot fomation). This model used 30 chemical species and-lOS reactions. 
The results were computed Using an adaptive *mesh placement strategy, 
and they resolve detailed structure within the flame. Note that the flame 
is very *thin''- its thickness is on the order of one millimeter, while com- 
bustion chamber dimensions are on the order of tens of centimeters. 

An interesting aside is to note that these two approaches to combus- 
tion modeling match corresponding approaches to experimental investiga- 
tion. That is to say, it is usually not possible to measure or pompute minor 
species concentrations in complex turbulent flows, where^^ it is possible 
' to cfp so in laboratory flames. 




FIGURE 2^5 Species mole fraction profiles showing the internal structure 
"of ai^ atmospheric-pressure stoichiometrfc^acetylene-air fiat^e. (From X A* 
Miller and'R* J. Kee, Sandia National Laboratories.) ^ 



The inherently disparate time scales of chemical reactions^ in other 
words their stiffness, contribute to the numerical difficult of solving com- 
bustion models (see Section 3.2). The inherently disparate time and length 
scales for fluid transport, heat transfer, and chemical reaction are respon- 
sible for the presence of steep fronts in the solutions. Also, there are many 
degrees of freedom in the ^stem of equations. The number of governing 
partial differential equations is large because a transport equation must be 
. included for each species involved.in the chemical reaction set. A system 
of 30 to 50 species, involving 100 or more chemical reactions is Epical even 
for fuels as simple as methane. Also, for practical combu^tors, the model 
must ultimately include complex three-dimensional geometries. Because 
complete models of real combustors are too large for present computers, 
ap important challenge is to simplify the models (including the physical 
assumptions) to a tractable level. In addition to these problems^ there are 
potential difficulties related to scaling. The ten^atures are on the order 
of 10^ K, while some species can have important efifects even Tvhen their 
mass fractions are as low as 10"^^. Moreover, before the computation^ 
. the peak ma^s fractions of the various species are usually .known only to 
within several orders of magnitude. " 

Many of th^ challenges of combustion modeling have^een met byJ;he 
^.numerical-analysis communi^; however, many more^^wait resolution. For 
example, for ^stems of ordinary differential equations wftgigdfica ^^d ! 
to treat Ihe stiffness that results from the compl^ chemical Icinetics by 
using stable implicit methods^ However, for systems of partial differential 
equations the application of these methods leaves many open questions 
about ho^ to treat the linear algebra and how to compute the error 
estimateSr Operator splitting methods are important in rendering the 
linear algebra tractable for large problems. Stiffness' also occurs in low- 
Mach-number flows due to very high veloci^, but low-amplitude; pressure 
wayes. Usually^ we do hot care. about the details of these, waves, but 
they can tinreasonably limit the size of the time step in explicit methods. 
Subcycling methods, rather than implicit methods^ are sometimes used, 
to alleviate this problem. Subcycling is a form of operator splitting (see 
Section 4.7) in which the invicid hydrodyntoiic equations are solved with 
small time steps, while the viscous parts of.the equations and the energy 
and mass transport equations are solved using, a much larger time step. 

Another particularly important topic in combustion models. conoe^s 
the adequate resolution of localized behavior such as flame fronts. Ohe 
line of research considers. adaptive meshing strategies in which a spatial 
mesh n^work is adjusted dynamicaliy so as to capture tlie local behavior 
accurately. Other work considers front-tracking methods, where the flame 
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is treated as a local discontmuity and conservatioii equations connect.botb 
sides of the front and predict its movement. 

Combustion models must account also for Rtud turbulence, an area 
where even the underlying physics is not well understood. Here, computa- 
tional' models such as the random vortex method are proving valuable 
in simulating turbulence effects, {Specially in the investigation pf large^ 
scale turbulent' eddies, the 3i>ca!led coherent structures. We classify thes^ 
methods as **problem-dependent methods'' because the physics and the 
numerical model depend heavily on each other (see Section 4.9). Recently, 
4he-randctn v o tt<jA mt*tlioU liM bee n combined with a flame propagation 
algorithm so that combustion events ean be modeled. Figure 2.6 is a se- 
quence of computer plots that shows the vortex velocity fields and flame 
fronts as computed from a model of turbulent combustion behind a step. 



2.2.1 Asymptotic Analysis 

We mention here that applied analysis (in concert with computation) is 
valuable for many problems in combustion. Asymptotic /methods, for 
example, can actually take advantage of phenomena such as steep fronts. 
They can thus be used to provide insights and to suggest approximations 
that help to reduce the complexity of the numerical models and thus render 
them more tractable. An important aspect of the asymptotic analysis is 
the possibility of considering the dynamical stability of flames. Such work 
may lead to fundamental^ understanding of such phenomena as the onset 
of turbulence. ^ 

Asymptotic methods ei^tqit the fact that the overall activation energy 
of the chemical reaction b typically laxge^ a well-known consequence of 
which is that flame fronts are very thin. That is, chemical re^tion is 
only important in a thin region where the temperature first approaches 
its burned value. On the nnburned side of this r^gion^ chemical r^tipn 
is negligible because the temperature is too low, and on the burned side, 
it is negligible because the reaction has essentiaUy gone to completion^ 
depleting the unbumed fuel. Mathematically spe^ng, this thin cbemi* 
,cally reactive region may be thought of as an internal boundary layer,, 
separating the unburned and burned gases. In the Uinit of a^mptotically 
larg^activation energy, the boundary layer ia infinitely thin, and we may 
use asymptotic matching principles to conne^t^ or match, the solutions in-^ 
. side and outside the boundary layen The result is a flame sheet model m 
which the solutions on either side of the^^eet are connected by nonlinear 
jump conditions that depend on local conditions at the/ront. Though the 
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resultibg problem is still nonlmear^ it.represents a significant simplification 
over the original problem with Anhenius kinetics and is in fact equivalent 
to the original problem in the asymptotic limit of infinitely large activation 
energy. 

The fint studies using large activation energy asymptotics began to 
appear about a dozen years ago and vr^re noteworthy for providing for the 
fint time analytical fkmne speed formulas for a steadily propagating planar 
flame^ The most significant results to come out of this approach to date are 
the predictions of instability and bifurcation phenomena, which describe 
celtutaf and pulsating modes of fiame propagation. These more complex 
solutions branch, or bifurcate, from the basic solution (in this case, a 




FIGURE 2.6 A sequential series of computer plots displaying velocity 
fields and flame fronts in turbulent combustion behind a step ^ inlet 
Reynolds number of 10000. The combustion is^stabilized by a recirculation 
region behind the step. The unbumed gas is a mixture of propane and air 
at an equivalence ratio of 1/2. (FVom A* P. Ghoniem, A* J. Chorin, and 
A. K. Oppenheim, Numerical modelling of turbulent fiow in a combustion 
tunnel, Phil TVarw. R. Soc- London! Ser. A 304, 303^32$ (1982).j 
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steady, f^btb^ flame) as one or more parameters exceed critical values for 
instability of the basic solution. Their significance lies partly in the fact 
that they represent intermediate modes of propagation of the Same in its 
evolution from a laminar to turbulent form of propagatiAi. Experimental 
observations of these transitional modes in their pure form have recently 
been reported, and the theory thus identifies critical parameter thresholds 
that separate me form of Same propagation from another. A nmnber 
of biftircation analyses using nonlinear perturbation methods have been 
successful in characterizing these nonsteady and/or nonplanar flames for 
parameter values near the instability threshold. However, numerical ap- 
proaches may be used to great advantage in describing these transitional 
solutions when parameters far exceed their critical values for instability of 
the basic solution* 



2.3 PLASMA PHYSICS 

Controlled nuclear fusion provides a possible, long-range energy source. 
As currently conceived, a fusion device must contain a deuterium and 
tritium plasma for a sufficiently long time for net energy production* The 
confinement may be effected either by magnetic fields or by simple iner- 
tia effects Inertial confinement of a pellet, coupled with laser or particle 
beam heating, is a distinct possibility, although the bulk oX^he world 
fusion program centers on magnetic confinement of a plasma. The un- 
derlying problems in magnetic confinement are the determination of the 
equUibrium, stability, transport, and heating properties of plasmas under 
realistic conditions. Numerical modeling and computation represent major 
tools in this study. In recognition of their importance, We Magnetic Fusion 
Energy Computing Center has been established at Lawrence Livermore 
National Laboratory in Livermore, California. This Center is the third 
largest scientific computing center in the United States, and its sole fimo 
tion i^ to provide computing capability to the scientists and engineer* in 
the U.S. fusion conmiunity. Fusion plasma physics spans a diverse coUe^ 
tion of fields, with significant Efforts occurring in the fields of Hamtlbooian 
particle. dynamics, statistical mechanics, kinetic theories, and dissipative 
and nondissipative single and nijiitifiuid models. The complexity of the 
problems involved mafces computation an integral part of the research and 
development program for fusion. \ '-^^ 

ControUed^fusion confinement experiments have indicated that, even 
in grossly stable configurations, fluctuations may play, to important role ia 
determining energy and particle transport. Thus, an unjderstanding of the 
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nonlinear behavior of various plasma models 19 fundamental to describing 
9uch behavior. For both ideal magnotoh}^drodynamics, a flui^ jnod et, and 
the Vlasov-Maxwell system/ a kinetic moddj linearized equnEhs ha\re 
been studied extensively. These ana^ses are appropriate for oescribing 
small^amplitude deviations, from a quiescent equilibritim but omit the 
effect? of mode-mode coupling and the onset, properties, and evolution 
of turbulence. The development of a self-consistent model of plasma 
equilibrium with fluctuations, stochastic particle and magnetic-field-line 
b^avior, and resulting transport will continue to be an impcfta 
area. 

Plasma fusion applications present many problems in which the equa* 
tions of Hamiltoniaft, dynamics appear. These systems describe single- 
particle phase^space trajectories in the Vlasav-Masnvell theory, of colli- 
sionless plasmas. In addition, the trajectory of a magnetic field tine in 
a toroidal system is described by a nonlinear Hamiltonian. Thus> the 
question of the existg^ce and construction of aiiiabatic invariants, ex- 
plicitly time-dependent or notj is ofnindamental importance for these iy&- 
tems. Significant progress in this area might be made by combining ideas 
from mod^ topological dynamics, numerical simulation^ and perturba^ 
tion analysis. For example, if magnetic field lines are ergodi(^ throughout 
a volume rather than lying on closed invariant surfaces, as^ia^^ven by 
the Kolmorogoff-Amold-Mpser thjeoryj then, owing to electron streaming^ 
the thermal conductivity within this volume will be very fast. Tbus^ such 
ouestions,as when ergedicity arises and what the properties pfHamittotiian 
djiTtamics ar^ under^ergodic circumstances impact strongly m these ap- 
plications. ' . 

Since computational modeling of the full three-dimensional plasma 
equations is out of the question, current work utilizes a range of different 
compromises, simplifying to di£ferent degr^s* ather the physics or the 
geometry to obtain a number of computationally tractable problems, each 
of which illuminates In a distinct fashion^ a different aspect of the full 
plasma problem. In spite of these diiBcultieSj the computational approach 
is a msjoT route to progress on th^ problem of controlled fusion. The 
reason is simple. Experiments are expensive^d must be supplemejited to 
the i^fiaximum extent possible by theory. The ^heory Is highly complex and 
nonhnear and is obtained by a combination of numerical experiment^ and 
physical intuition. Clearly improved computational methods are one of the . 
methods through which progress is achieved in this area* Development of 
more powerful computers will also be. requured, as it is hard to believe that , 
smart algorithms will by themselves suffice; ' 
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2,4 PARTICLE PHYSICS 

The problem is to find the equations that describe the elementary particles 
of subnuclear physics. This problem is among the most difficult to have 
been considered seriously by science in this century. In recent years there 
have been significant gains in our understanding of the mathematical 
structure of the equations of quantum field theory. An analysis of the 
mathematical existence question for the quantum equation 

shows that in space-time dimension d< if the theory exists, where^ in 
d> 5^ there is no such theory (as ivould be given by standard methods) 
except for the trivial case X = 0, The method- of proof sugg^ts nonex- 
istence for the physical case d = 4 as well and offers ccwfirmation from 
the side of exact matheanaticAl analysts of ideM"a5vahcea~By"tTrebretical' 
physics. This mathematical confirmation is much more compelling than 
any confirmation yet ofi^red from experimental physics. 

The modd is ^^near^' (o the Sine-Gordon equation, which is com- 
pletely integraUe^ Among recent results are the discoveries of a large 
class of stable, localized (in space), periodic (in time) solutions (analogs 
of the Sine-Gordon breather) and of the existence of breather formation 
resonances in soliton/ antisoliton scattering. There are many relataf chal- 
lenging problems in nonlinear iQath^atiCs and dynamical systems theory. 

The problem with the equation in c( = 4 dimensions is one of several 
reasons for- considering in its place quantum gauge^fields. A standard 
approach to the quantization problem reduces it to the existence question 
for a singular non-Gaussian functional integral over an infinite^dimensional 
space. In the ease of^gauge fields^he gauge potentials form an infinite^ 
dimensional afSne space A on which the group § of gauge equivalence 
* acls< The fi^ctional integral is only defined over tlie infinite-dimensional 
manifold 

^ " M=A/g 

This manifold is not Qat The integration over M has been shown by exact 
mathematical analysis to require introduction of coordinate patches. Only 
locally, within a single cppfdmate patch can M=A/^he regarded as an 
open sutset of Euclidean (Hilbert) space. 

One of th^ methods proposed for understanding integration ove( ^ is 
to understand the critical points A6 M. This question leads a study of 
classical solutions of the Yi^g*Mills equation and a reduced form of this 
equation, caUed the self-dual Yang'Mills equation. Here we are looking. at 
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a specific nonlinear etltptic equation in d — 4 dimensions.' A remarkable 
analysis has led to a complete classification of its solutions^ using methods 
ofgeometryr topology, and al|^bra (fiber bundles^ the index theorem, and 
'filgebraic varieties). ^ 

The developn^ent of ^liimierical methods for the 3tu(iy of quantum 
fields is in its infancy, mainly because thcproblem is to a large extent out 
of range of present-day computing machines. Among the methods used 
on this problem^ we mention Monte Cai!# integration over a space-time' 
discretized version of X = AfQ, The discretized problem is called a lattice 
gauge fields and even after discretization^ the 'dimension of X is too high 
to allow evaluation of integrals by direct quadrature. Asecond approach is 
to generate series coefficients from cluster expansions of statistical physics 
and to reconstruct the desired function space integral from a Pad^ analysis 
of t&& coefBcients. Other attempts-have been based on finite elements and 
on the renormali^ation group. It may be that a collaboration of numerical 
analysts with mathematical or theoretical physicists oikthis problem would 
■"hvbeneftcifil:^^^^ — ^ ■-■ = — r 



2.5 CON]!)ENSED^MATTER PHYSICS. 



2.-5.1 Statistical Phyaics 

In this subject, the equation of state, transport coefficients (pressure, 
viscosity, aad thermal conductivity, for e^cample), ^nd other n^acroscopic 
properties, of matter are related to and derived from intermplecular forces. 
The area has diverse and important appUcations ^ ranging from metallurgy 
to pQlymer chemistry to semiconductors and is an active area of research 
from the points.of view of theory, numerical methods^ and experiment. 
Here we limit ourselves to the two mathematical aspects of this sabjeet: 
mathematical 'theory and humericat^methods. 

The equations of statistical physics^volve a large or infinite number 
of degrees of freedom, and so^e mathematical theory^of use here is 
analysis over idSnite^imensional spaces. Almost the same mathematical 
structure arises as in the study of quantum fields, the relationship between 
the two being th^t a quantum field js a continuum limit of a statistical 
physical (crystal^ lattice) model. ■ — 

^ In a small cl^kss of models, including the two-dimensional Ising model, 
exact 8olul(ons are known.^ larger but still restricted class of models has 
been analyzed mathematically with, respect to qualitative behavior. One 
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FIGURE 2.7 Monte Carlo simulation of a kinetic Ising model compared 
with x-ray and neutron defraction measurements on an alloy "of 60 percen^^ 
gold and 40 percent platinum that has been heated, and then quenched to 
60 percent of its critital temperature. The abcissa is a reduced momentum 
transfer^ and the ordinate is a reduced scattering intensity. (FYom M 9. 
Kalos, New York University .) ■ ' 
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issue recently ^dressed in this work was the stability of sulfaces^ as is 
relevant to the problem of crystal growth. Another is the effect of aperiodic 
or random crystal structure on bulk material properties. A striking recent 
development was the mathematiciil demonstration of the dipole finding 
(Koaterlit^^-.TEouless) phase transition. . - * 

The numerical methods for statistical physics are. basically of three 
^es. In molecular-dynamics calculations^ one takes a large number of 
particles and follows their motion by inflation of the ordinary difibential 
equation ■ . , . " - , 

defined by the intennolecular forces F^. This. method is,"exact^ in its. 
treatment of intennolecular forces, to the extent that they are known 
and that quantum effects can be neglected,, but it is approximate in its 
treatment of statistics, sinc^ there are computeP'dictated limits on the 
number iV of particles that can, be included. ^ ^ ' - ^ 

^ - "'"^^ 



Next iTO mention the method ^ series expansions, applied to the 
calculation of the equilibrium distribution and the partition 2y which 
typically has the form - a 

(ic;r=e-M^)n^ \ ' ^ ^ 

aod ' * 

where r is the intermolecular potential .^nergy* Then thermodyiffimic 
functions su^h as prepare emerge as derivatives or Ic^arithmic derivatives 
of Z with respect tp the parameters (e.g., in V). If the potential energy 

. ? ' ' . ' 

is a sum of pair potentials V then the identity 



= l + (e-^-l) 
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FIGURE 2.8 Energy of liquid helium as a function of density. The broken 
line and points with errors show results of the Monte Carlo^quadrature of 
the many-body nonrelativistic Schrodtnger equation. Th^ solid line is a lit 
to the experimental data*,{IiVom M. H. Kalos, NeW York University.) 
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substituted in the definition of du) above yields an e^cpansion that is rapidly 
convergent^ where ^ 
. . |e"^-l|<SCl 

i.e., where U^O. ^his is the region of noninteraction, and so series ex- 
pansions'are especially useful to give weak corrections, for example, real 
gas corrections to an ideal gas. The expansions appear to converge up 
to critical points and with considerable work have' been used to extract 
information on the equation of state in that region. 

. Monte Carlo methods are used in the quadrature of very large- 
dimensional spaces, such as the determination of (Lr above, as well as in 
the direct simulation of stochastic systems. These methoda have become 
an experimental tool of mathematical physics. In studying the qualita^ 
tive an d^ quantitative behavior of a highly idealized model such as the 
Ising modeljp equilibrium or very far from equilibrium^ or of lattice gauge 
theories, cohtinuum and lattice models ofj)olyniers, and atomic models of 
quantum liquids and solids, one can carry out numerical studies in ^Vhich 
the high-dimensional (i.e., many-body) character of the' problem is not 
distorted. 

In appUcation to the ^calculation of the equilibrium distribution (Lr, 
the essence of the Monte Carlo method is as follows: Starting from an ar- 
bitrary point in a given ensem^, one n^pdifies a single particle position 
"at random^ but usually so as to lower the potential energy V, occasionally 
so as td raise V. After enough such elementary dteps, convergence to the 
distribution du) is obtained. As describe here, the method is extremely 
simple, and complications arise from the necessity to obtain convergence 
in a reasonable time for realistic problems. 

- Perhaps the most significant success has been'in the microscopic 
theory of classical Buids, where Monte Carlo mojdeling has provided the 
"experimental** basis for the accurate expansion int |e"^ — 11 mentioned 
above. Another conspicuous success is the extraction of critical exponents 
by joining ideas of the renotmalisation group to Monte Carlo simulation. 
Figure 2.7 showsMonte Carlo simulation of a kinetic Ising model compared 
with scattering measurements on aieal alley. 

Figure 2.8;shows how Monte Carlo methods applied'to the many- 
body nonrelativistic Schrodinger equation give a quantitative account of 
the energy of ^real) liquid helium as a function of density. 
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2,6 GEOPHYSICAL ^PUCATIONS 



Out of a wide rang? of geophysical applications yte limit ourselves here to 
threci modeling^ problems iii connection with the discovery and production 
of petroleum Mathematically the discovery process is an inverse problem 
(see Section 3,6)-that of constructing geologic maps using seismic signals, 
which in turn are generated by vibrations'oj; explosions at the surface of 
the Earth, As the signals penetrate the Earth, they cross layers of differing 
density! This causes reflected signals to be returned to the surface, where 
they are then recorded and ana^zed. The raw data are very noisy^ owing 
to irrelevant near-surface density fluctuations. The noise is removed by 
averaging signals from neighboring receptors or from neighboring source 
locations. Then multipje^flections mustiBe subtracted^ and a^Qpen- 
satton (called normal moveout) must be introduced for effects ottionver^ 
tical signal propagation. Reflection from nonhorizontal. lasers ^Herates 
complicating shear waves as well as pressure waves. The subtraction of 
multiple reflection signals can be based on Foiuier analysis an a half^space 
and a Wiener-Hopf facforization. 

The correction for nonverticuT jaignals can^J>e based on a reduced 
Helmholtz equation, (-A + fe^ + v)u 0, but in brder to focus on flig*^ 
nals moving in only one direction (either up or down)» one takes first' a 
square ipot and then a ppwer series expansion" of the square root. This 
process is, known ^ the pmbolic approidipation ^d leads to the famihar 
Schrodinger equation. Alternatively, ^he analysis can be based on rsy 
tracing and the solution of ordinary differential equations. The numerical 
problem is to implement these steps efficiently, in view of the large amount 
of data to be analyzed. 

The production problem is to describe the flow of oi]» water, chemicals^ 
and/or heat in a porous sandstone layer, ' The equations to be solved 
are a coupled system of nonlinear parabolic equations. Generally- the 
equations describe mass ^nservation of individual species^ and typically 
some are nearly elliptic whereas others are nearly hyperbolic. In the 
hyperbolic equations, coherent shock and rarefaction waves describing oil 
banks^ absorption fronts^ and flame fronts may form (see Sections 3,2 and 
3,3), Depending on a dimensiohless numbet known as the mobility ratio^ 
heterogeneit^^esj and geometrical effects (convergence versus divergence)^ 
the fronts may bekome unstable with respect to the formation of flngera 
(see Section 2,i). 

Critical issues in this problem are the eflicien^Tsolution of large sparse 
linear systems arising from space and time discretizations and the con^ 
trol of numerical instabilities in the solution methods. The^hysical in* 
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FIGURE 2.9 Successive time steps in the mofvement of an oil-water in- 
terface toward a ptoducing well The view is a vertical cross section, 
with the well locat^ on the left boundary of the computational region. 
The^stabitity is caused by three factors; by a heterogeneous conducting 
channel near the bottom of the reservoir, by the converging cylindrical 
flow^ttem near the producing well, and by the better flow properties of 
the less viscous displacing fluid (water). fVom J, Glimm, B. Lindquist^ 
0. McBryan, and S. Yaniv, Statistical fluid dynamics E: the influence of 
geometry on surte^^ iustobihtiest in FronUers in Applied Ma^ Vol. 
,1, R. Ewing, ed.,]gIAM, Philadelphia, Pa. (to be published).) 

stabilities mentioned above are controlled ""by the use of heat (to ihake 
the oil flow more easily), or the use of a heavier pushing fluid (polymer- 
thickened water in place of water or water and CO3 in*place of CO2). , 

The third modeling problem we discuss is the process of Uftang th^ 
oil to the surface of t'h^ Earth. Because of the sudden drop in pressure, 
considerable gas will typically come out of solution, and the resulting two- 
phase mixture is in the droplet or niist regime mentioned in Section 2.1 . 
In contrast to the case of reservoir flow, even the basic equations for two- 
phase flow in a pipe are;not well understood. 
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2.7 METEOROLOGY 



■ Accurate forecasting of the movement of large-scale weather patterns is 
clearly an important problem. Here we mean the trackin^of highs and 
lows f>n the order of 1000 kilometers in diameter. In addition, we are 
also concerned with the prediction of cloud patterns, precipitations, teip- 
»perature, and other elements of the weather in as much detail as possible."^ 
Global predictions of the atmospheric flow are now made routinely at as 
many as 12 levels in the troposphere and stratosphere.' . 

By studying dynamic models of the atmosphere, estimates have been 
made of^^the sensitivity of th| atmosphere to small perturbations, often 
4 quantified in terms qT the error doubling time. If *two initial states of the 
atmosphere differ by random variations of 1^, then it is found that ^e 
resultant states will diffi^ by 2*^ in about 3 days and by about 4*^ in about 
6 days. Hence, meteorologists have said that the weather is unpredictable 
for periods of more_than 2 to 3 weeks> But, the forecasters are still far 
from being able to predict, the weather accuratSy for even 1 weeE. The 
best results are currently obtained by the European Centre for Medium 
Range Weather Forecasting, whose 4-5 day forecasts are currently superior 
to those of the National Weather Service. 

The accuracy of short-rafige predictions is limited by four somewhat 
distmct factors: ^ 



1. Initializations: (a) Th^ accuracy and completeness of obserfational 
data and (b) the compatibility of these data with the mathematical model, 
which is, of course, a simplified repre^ntation of the atmosphere. 

2. The limited number of grid points or spectral components used 
in the truncated, i.e., numerical, model as a result of limited computer 
power. ^ ^ 

3. Missing^or severely truncated physical processes such ^ cloud 
dynamics, prog&r representation of turbulence, and radiation fields, lead- 
ing to inade<{(fac]^ of the mathematical model describing the atmosphere. 

4. The inherent finite limits of predictability of certain types of non- 
linear dynamical ^stems such as the at»osj>here. (Compare ^he previous , 
paragraph.) * - 

^ \ % 

The study of problems related to the reduction of the influencen)f 

these factors is at the heart of^most research that is directly related to the 

improvement otweather prediction. - 

The sparsity and inaccuracy of observational data 11(a)] are being 
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partiafly ov^come bjrthe use of weather-observing satellites. Concerning 
1(b), considerable research is currently directed toward the analysis and 
filtering of the observed dat^a^so a$ to mak^hem compatible with numeri- 
cal models. Items 2'and 3.are closely related==!fiinc€, in a general way, the 
lack of a coDCkplete representatipn of all physical processes is related to^ the 
incomplete numerical resolution of the mathematical model. Such physi- 
cal processes include cloud dynamics (and the associated cloud physics), 
boundary-layer dynamics (including detailed features of the terrain and 
local heat, moisture, and dust sources), and detailed radioactive processes. 
In item 4 a new concept, "finite Umits to predictability" has emerged 
froM^the work of E.'N. Lorenz, namely the idea that the fitmosph^Te.and 
perhaps certain other dynamical systems majr have finite limits of pr^c- 
t^bility regardl^s of the accuracy and detail of their initial conditions 
and the accuracy with which they are computed. Lorenz concluded from 
a simplified model^that hydrodynamical ^sterns with an energy sp^trum 
(in the customa^ sense) having- a power*lawexfH>nent greater than —3 
(e.g., -5/3) should be lihpredictable in, detail ^SUt a finite time. He sug^ 
gests a time of about iS^days for the atmosphere, whose large*scale energy 
sjpectrum seems to haVe an exponent of about? -2.8, with smaller scales 
(i!e.,TbaB than say, 20^1an) Having the -5/3 power law associated with 
homog^eous turbulence. Thi| questicm of predictabUit^ is clearly related 
to that of the stabiUtyjof hydrodynamical ^^tems and to modem th«>rie8 
of bifurcatibn and chapB (see Section 3.5). ' 

A large <^omponehi of atmospheric research is concemed with topics 
that have tittle direct application to forecasting. Their motivation may 
lie in basic scientific inquiry or in appUcations to other disciplines, e.g.,, 
radio-wave transmission for satelUte communications. They may involve 
the study of fundamental hydrodynamic phenomena or other physical 
processes several steps removed from specific application. In all of these 
studies we find extensive use^of analytical and numerical modeling of 
s^hysi^ Bystemrand fiffl ce^. .. In. Edition, a yi^de varie'^ of statistical 
approaches are used in all are^ oT^^tmospheric research. ^ ^ 

Let us mention a few mpre 'ai^e^^ studjf. In addition to the praj^ 
tion of large-scale weatfier patterns.^computers.ai^ebm used in.Umited^ 
regions to study the devefopn^ent arid tooifem^tp^ 
bancea,.|3^h as hmriqah^, thunderstorms, and tomlSioesTHei^dpptive 
griti^methods are used\(see Section 4.2). It is hoped, fpj^examp^that a 
better uiiderstanding^of the detafis of these destructi^^phenomena will 
lead to thepossibili^^of altering their course and deydopinent. ^ ^ 

<Sn the longer tiihe scale, where seasonal forecasts ofaverage rainfall 
and temperature ar^-made over large regions of the globe^ jthei^ ^ at^ 
present several potehtiaily useful, but .unpiroven ideas*' Here, however/ 



some r\evi mathematical and physical insight is needed in order to develop 
a satisfa(;tory "average" system of equations^ This is an importantrx^pen 
problem. - * 

rOn the^ttU longer time scale of decad^and beyopd^ d^rmining 
the effect of our mechanized civilisation on the environment is a bask 
problem* Significant progress has already^ been made, through the use 
of simple climate models in conjunction with paleontology, astronomy, 
gpology, and volcanology, in efforts to understand the factors that influence 
irregular alternations of ice ages and interglacial periods. These statistical- 
. dynamical climate modeb will be essential for predicting the long-term 
effects of the observed increase in carbon dioxide and similar problems 
that will almost certainly arise as our industrial civilization expands. 



**2,8 ASTROPHYSICAL APPLICATIONS ^ - " 

Astrophysical studies pose a wide range of problems for mathemati- 
cal and computational analysis. The underlying physical theories in- 
clude hydrodynamics^ nj^etohydrodynamics and plasma theory, radia- 
tive transport theory, atomic and nuclear physics, and general relativity 
theory, and hence various comments made in <^er applications sections 
^ of this report apply this area- ^Hfe. 

In recent years considerable progress bas^^^^^e in theoretical and^ 
numerical studies Qf a number of astrophysicaltopics^as, for instance, tUxe 
study of stellar interiors, the forination of stars and galaxies,, the spiral 
structure of galaxies, the physics of supemovae and the evolution of su- 
pemovae remnants, the formation of the solar system, and the behavior of 
binary star systems, to mention Just'a^ew. The computational approaches 
may range from relatively simple simulations to the numerical solution of 
complex systems of partial differential' equations. - 

As an example of the 6rst type, simulations have been applied to 
the study of star formation in a galaxy. It may be assumed that vthen a 
massive star becomes a supernova the shock mve emanating from it can 
compress the surrounding ^terstellar gas creating nejy^stars. If at least 
one new star is also a massive^star the phenomenon can repeat, leading 
to a chain reaction in the creation of stars. This is equivalent to a direct 
percolation problem, and, as is typicalfor such problems, phase transitions 
ate involved. This is a nonlinear problem with a complicated structure in 
spaqe and time. Accordingly, analytic techniques are difficult, whereas 
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computational simulations are relatively straightforward and lead to the 
development of realistic model galaxies^ ' 

M a second type of exaniple vre mention a problem that is amenable 
to considerable mathematical analysis as well as to computational attack. 
This is the question of the spiral stnictxire of galaxies. The so-called 
density wave theory approaches this question as a dynamical ^problem in 
the form of the gravitational instability of a galactic dish wi th, respect, to 
spiral modes. There are three basic approaches in tiie calculation of spiral 
^ structures of galaxies on the basis of density wave theory that may be 
charg^rized loosely as the stellar model, particle model, and Suid model. 

The' stellar model represents the classical ^proach in the study of 
stellar systems. In brief^ a galaxy has a stellar component and a gaseous 
component. The latter usually has asuffidently small mass to be negligible 
in first approximation. The baac equations .governing the behavior of the 
stellar component are the Boltzmann equation, usually in collisionless form 
when close encounters between stars are omitted^ and the Poisson equation 
for the gravitational potential. In a sense, these equations are much simpler 
than the Vlasov equations of plasma physics^ but their numerical solution 
nevertheless poses chaUenging questions* 

In the paiiicle models the stellar component is considered as a very 
large but finite number of particles. As in our first example, the computar 
tional approach then assumes the form of a simulation process. In brief, the . 
motions of the individual stars are followed, and their gravitational fi^d 
'is calculated by a self*consistent evaluation of the field. This technique 
has provided valuable qualitative information But the number of stars in 
these models is usually of the order of 10^^-^^ and hence^ by necessity^ 
pumerical simulations are extremely limited in accuracy and provide only 
few quantitative data^fos^^cified galaxy models. 

In the third approach^ the 3tella!r system is considered to be a con-^ 
tinuous medium. In this setting one may ttien study the characteristics 
of wave patterns over the galactic disk and the n^ dependence on the mass 
distrib.ution. In particular, the spiral structure of galaxies may be ex* 
plained in terms ot spiral wav^ patterns of some kind. Such^ stationary 
wave patterns over a field of difierential rotation are to be expected; in 
fact, l^rodynamic waves over ^hear fl<>ws are kno^yn exist for some 
time especially in the form of seif-excited modes. The general form of ^ 
the resulting theory of spiral galacticjtnictures requires elaborate com^^ 
putationai t^lmiques. M the same^me^^, asymptotic approaches have 
■provided, analytic results that haiw Jed to a better understanding .oFUie 
dynamical mechanisms de^ite their limitations m accural^ and in^thiet 
types .of galaxies ccrored, , - . , ^ 

As noted earlier, these are only some examples of the many diSerent 



types of problema in astropbyaics that relate closely with applied and com* 
putdtional matbematics. Many of these problems involve extremely wide 
rang^^ length and time scales. For example^ calculations of^stellar evolu- 
tions have to range over billions of years, while the dynamics^ say^ of the 
supernovae phase takes placf^ vdthm iniUisec<xids. Similarly, the size of a ' 
-neutron star-diffeK-by substantial orders of magnitude from the size of its 
corresponding relevant gravity tteld. i hus^ m piamcuiar, the mAthem&tlc&l 
and computational difficulties^ discussed in Section 3.2/are especially ap- 
plicdbie here. Moreover, the wide range of the underlying physical theories 
leads to many substantially diSeient types of mathematical models, which 
in turn require very difierent computational techniques. 



2.9 STRUCTURAL MECHANICS 

During the past two decades, the use of computers has transformed 
large parts of solid mechanics into practic^ toob, for ^ multitude of 
techndogica] developments. . Sopbisticat&i ,cc^pu4ieLU0^ so^ is 
employed^hroughout the nation's industries and resei^ch laboratories in 
the analysis and design of structures and mecbanii^^quipinent There 
is a strong, interacti6n-between applied mathematics apd solid mechanics. 
Matbematicat analysis has provided insight into model formulation and 
the development of powerM numei^c^ methods; and^ vice versa^ novel 
engineering approaches have led to new research areas in applied mathe- 
matics. 

In the case.df ^linear problems there, exists now a relatively broad 
experience in computing solutions for a. range of problems concerning 
the behavior of solid bodies subjected to Specified loads. In general^ 
the computed results are considered reasonab^ reliahle, and they have 
been corroborated over a period of time by observation and practical 
experimen^tion. There appears^ however, to be a growing need for the 
develppment of computable error estimates that can provide a realistic 
cfieck on the solution accuracy. Such a posfcenon bounds have been shown 
to be feasible; but their applications to large, realistic problems in solid 
mechanics still requires considerable research and software development. 
The general availability of economical a pos^non^stima^^vt^^ make 
po^tble the consistent use of adaptive m^'^refinenient^£e<foniques tHat 
would reduce the cost of data preparation by the users and make it possible 
to generate near^optimal solutions for ji given ainount of computational^ 
expense. ^ 

my . # 
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The general trend in computational solid mechanics today is toward 
extending the computational methodology to nonlinear pp)blem3i^ Soui:ce3„ 
<3t nonlinearity in structurial problems are (1) geometnc|nonlmeiaity due 
to nonlinear strain-displacen^t relations, (2) material ponlinearity due 
to nonlinear c<mstituti;ve equ^ions, (3) force noolineari^^&due to nonlinear 
stress boundary'conditionaj. aiid (4) kinemati c constraiy t j^^ <\\\^ 
to'honlinear di$^er^ent ^imdary conditions. The sofi^ of nonlinearity 
affects ihe Torm^r the^resulUiig nonlinear equations and^^ence^ influences 
the effectiv^e$5^th4 soltltion techniques. 

The numerical ^alyeis of all of these nonlineai^ problems is not yet 
at a satisfactory stage. Many computer programs Tor such problems exist, 
but the mathematical pms Tor most or the methods used is insufficiently 
understood, and liier^s Ihtle known about t1^accura<^ or the computed 
results. Nforeover* in tjie case of nonjinear probbems, few numerical com-, 
putations can be ^up^emented mth siiScjent experimental experience. 

The situationMs best understood in the case oT finite elasticity. 
Even there th^ mathematical theory of the underlying equationsas incom- 
plete, and the approximation th^ry for these Equations is g merally based 
on various simplifyin^^lssumptions that may or may not :>e valid for a 
particular problem; \> " ^ ^ ' ■ ^ 

^,The state of l^he^ in elasto-dynamic problems is i: even wrse 
shape. It is^knovm th^^ultiple solutions may exist and sb«ck \vaves can 
develop. Moreover, n^fi^^ solutions are necessarily physically i^levant. 
The.cpiestions of how to model Such phenomena numerically and how to 
detennine the pbysicaUy realistic solutions are as yet largely open. When 
it comes to problenis in finite plasticity even less is known. Although there 
has been much progress in this area^ during recent years^ no satisfactory 
and complete mathematical model is available as yet Especially, there 
'are profound mathematical and computational difficulties in modeling 
please changes,, viscous effects, cracks and sing^arities/ the growth of 
cracks under dynai^iypadi^gf and the identification and implementation 
of pl^sically reasonabWcojrstitutive equations describing th^ materials. 

^aO NONDESTRUCTWE TESTING AND TOMOGRAPHIC ) 
-REQbNSTRUCTIQIl?"*^; - r 

Ther^ is a pervasive need in technology to evaluate quantitatively the 
integn^ and the rem|^&ig reliable lifetime of components and structures, 
e.g^,^rom bridge girdeff to higK-perfo^aiice ceramic turbine disks. In 
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the past decade coDsiderable progress has been made. This developing 
techDology is called noD destructive evaluation (NDE) to distinguish it 
from older noDquantitative Dondestnictive testing practices, NDE presents 
challenges to the applied mathematical sciences on mmy levels, ,A few of 
the more'important and topical problem areas are mentioned here, 

In NDE applieat i on s L a component is ofttat -s ubjectcd to som e .wt t uf 

penetrating radiation with the aim of deducing information about its in* 
temal state from a measurement of the radfation field external to the part. 
Examples include the use of ultrasonic radiation, x rs^, and neutrons. 
Because of its Sexibility^ relative cost, and safety, ultrasonic methods are 
often used in NDE applications. The deduction of information about fiaws 
from the incident and scattered ultrasonic fields relies on the solution or 
approximate solution" of the inverse scattering problem for elastic ^aves. 
In some regimes it is possible to develop and adapt i|na^ng techniques 
to the ultrasonic setting. In either case information dboiit defects must 
be obtained from band-Umited nois/data. A Hmdatnehtal limitation to 

''our current "ibility to utiUze ultrasonic techniques bfo^^ is limited 
understanding of the elastic inverse, scattering problem in either the fr^ 
quency or the time' domain, . ^ 

Jn otbe^important NDE applications^ efectrical currents are induced 
in a material. These currents produce gelds that vary depending on 
whether a defect, e^g.^ a crack, is present in the material. The utilisation 
of these so-called eddy-current^tediniques depends on the ability to infer 
information about the defects from the Treasured fields they produce. This 
again is an inverse problem that is imperfectly understood. 

Passive methods in NDS are also widely used* For example^ when a 
pressure vessel or aircraft is in service^ a crack once formed may grow and 
propagate. This will be accompanied by the release of acoustic energy. 
This energy can in turn be monitored at selected sites. The problem then 
beconies one of identifying and classifying the sources of these sound patr 
terns using acoustic emission studies, ^ increasingly important technique 
in NDE, In mathematical terms the^probl^ is the so-called inverse source 
problem, which is beset with the same 1y|>e of difficulties that the imrerse 
scattering problems pd^sses* Sparse, noisy data often taken at higfa^ 
noncptimal k)cations are the raw information from which source charac- 
teristics must.be deduced. 

The techniques df NDE have application in other areas, j^d much 
can be learned in other applications fields that is valuable to NDE, For 
example,, there is a^ck)se connection between fhe need for inverse scattering 
results in geophysics and NDE and acoustic ima^g . results inJ^E and 
biomedicat applications. An area of great success in medical applications, 
tomography, can provide? useful informaiion^in selected NDE applications. 
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The greatest success of tomography has been in medical applications 
such as the CAT (computer assisted tomography) scanner. Unlike ordi« 
nary x-ray technique^ which masks important features by superposing on 
a single picture infonnation on planes perpendicular to the beam, com- 
puterized ^c-ray tomography provides Pictures of indiyidufit thli> flli^^^a ' 
through the l^ody. Several hundred parallel x-ri^ pencil beams are 
projected in the plane of the slice^ and the attenuation of the x rays in 
each beam is measured separately and recorded. The procedure is repeated " 
for many different beam directions. An elaborate calculation then permits 
approximate reconstruction of the x-ray attenuation density as a function 
of positioii within the slice! 

The idealised mathemati<^'problem is the reconstmction of a func- 
tion of two variables from its integrals along lines. This problem^ as well 
as ^s'three-dimensional Version, was solved by Radon in 1917 and later 
re<Kscovered In various scittlngs such as probability theory (recovering a 
pfobabitity distribution ttora its marginal distributions) and astronomy . 
(determining the velocity distribution of stars from the distribution of , 
radial velocities in vanous directions}. Of course^ much work was needed 
to adapt the Hadon inversion formula to the' incomplete information avail*- 
. able in. practice. Various algorithms for the numerical inversion of this_ 
ill- posed proUem have been proposed, with the present trend favoring the 
S4>^caUed convolution algorithm on account ofits speed and accuracy. Each 
area of applications has, however, its own requirements and may need a^ 
modification of existing reconstruction algorithms or even a custom'-made . 
one. Some algebraic methods^ for instance^ can easily incorporate a priori 
. information about the object to be reconstructed. 

Keceiit advances in medical tomography include nuclear magnetic 
resonance (NMR) tomography and pqsUron emission tomography (PET). 
In NNIR, strong jnagnetic fields are used to affect the nuclear magnetic 
spin rate of hydrogen atoms. By varying the fields and their direction, the 
plane integrals of the density of hydrogen can be measured and the den- 
sity reconstructed' by an algorithm based on the above three-dimensional 
version of Radon*s theorem. The technique is now regarded as competi- 
tive with x-ray tomography for many purposes; and, of course, it is not 
ionizing. The advantage ^f^PET^ over CAT is that metaboUc processes 
can be followed. A^ compound such as glucose is made using carbon-11 
atptps, which einit positrons. Photons resultii^g from the annihilation of 
an emitted positron with an electron are detected by a bank of detectors 
that can record coincidences. Kecentty, algorithms based on probabilistic 
arguments have been proposed for the PET reconstruction problem. 




2.11 MATHEMATICAL MODELS IN THE BIOLOGICAL SCIENCES 



The biological maBifestations of the physical lawa of the universe pr^nt 
us with a rich variety of new phenomena that require the devlopment 
6r n^w m&th€th^Xic^ tdols and (:;6mpul8Ci^hai~^letho3s. We shall discuss 
Just a few examples of mathematical research in cardiovascular physiology 
and neurophysiology, with the knowledge that there are many other areas 
of biological sdences in whjl^ch mathematics and computing are fruitfully 
applied. 

Blood flow in the heart.obeys the incompressible Navier-Stokes equa- 
tions, which, in turn, are simply a statement of Newton's laws in 
differeotiai form (see Section 2.1). The distinctive biological character 
of the problem comes, however, from the moving boundaries that are in 
contact with the blood. These include the muscular heart walls and the 
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FIGURE 2.10 Computei^generated plots showing i^ne predicted opening 
movement of a^ curved pivoting-disk valve moun^ in the mitral (inflow) 
position of the/iSSft ventricle. The curvature^n^es the valve open more 
widely than a^stiaight valve piVote^ at the same point. It also helps- 
to prevent stagnation in the smaller opening of^the valve. .\From t>.,^. ^ 
McQueen and C. S. Peskin, Computer^assisted design of pivoting^disc^ 
Ijrosthetic mitral valves,,/ Tkorac- Cardiovc^c^ Surg, (in press).] 
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FIGURE 2J1 lYansmission and "reflection" of a nerv^ impulse at B Junc- 
tion where the diameter of the neuron sudde^y increases. The plots show 
computed voltage as a function of time at equall^r spaced positions. The 
Junction is at x = Xo, and the ratio of diameters ther^ is 2.5:1. Note the 
increase in propagation speed for x > Xq. A reflected wave iB set up .when 
the larger flber re^excites the small^ fiber after the refractory period of the 
smaller flber has elapsed. [From S. S. Goldstein and W. Rall> Changes of 
action potential shape and velocity for changing core conductor geometry^ 
Biophysical J, 14> 731 (1974).) ' . ^ 



elastic heart valve leaflets. The motions of these boundaries are not known^ 
in advance; they must be computed along with the^^tion of the fluid. 

These considerations have led to the develQpment^of a computational . 
model of the left heart that can be used in the cornputer-aided design 
of artificial heart valves. In Uiis model> the fluid ec^atlons are solved 
by finite^diSTerence methods on a regular^ square mesh (see Section 3.4). 
The boundaries are represented in Langrangian forrp as a ^collection of 
moving points. Coupling coefiicients between bounda^^markers and fluid 
mesh points are computed with the aid of an appromiation of t)irac's" 
6-fimction. This computer model (Fig- 2.10) has been used in the design 
of prosthetic heart valves to remedy problems of stagnation and blood 
dotting in the smaller opening of the valve. The i^bdd has also been 
helpflit in the study of disease processes, providing^ foT.example> a possible 
explanation for mitral valve prolapse^ ' ' 
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Just as the flow of blood in "the heart is ultimately governed by 
Newton's laws^, the conduction of electrical signals alon^ nerves is uK 
' timately governed by Maxwell's ^equations (Figs. ' 2.11 and 2.12). Here 
alsn ^ nature iffQvides a peculiar boundary caadition that leads to entirely 
new phenomena. In this case the boundary condition comes from certain 
voltage-dependent channels locate'^ in the nerve membrane. These intro- 
duce a noblinearity^ and the equations of nerve conduction take the form 
of nonlinear diffusion equations-tbe Hodgkin-Huxley equations. Without 
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FIGURE 2.12 Computed electrical sifn^s at .thfe,inlpui i&nd of a neuron* A 
^rief pulse of current is'applied at the periphery of a tree^ and the result- 
'mg voltages are computed {logarithmic 'scale) at the input (BI), at succe^- 
siye branch points (P, G?,, GGP), and.fmaliy at the cell pody (SOMA)* 
(Modified from J. Ruizel and W.'RaII> TVansient response in a dendritic 
neuron model for current ipjected at one branchy Biophysical J. li, 759 
(1974).] . ■ *' ' 
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the nonlinearity, signals introduced at one end of the nem w(Juld decay 
rapidly; because of the nonlinearity such aignals evolfe into a spectfipf 
wayeform tliat propagates at a constant speed withpu^ distortion. This* 

distance communication in the nervous system. 

' Currently^ there is a^viderahge of mathematical^ c^putational, an(| 
physiological research activity related to the Hodgkin-Huxl^ equations. 
Mathematically^ there is extensive research on the ba^ic theory of .non- 
linear diffusion equations, A {j^icularty fruitful approach h^e has been 
the use of piecewise Unear models that expose the basic structure of the^ 
equations. Singular perturbation meUfcds have also been useful because 
the equations exhibit a disparity of^pe scales. 

An important physiological ent^prise is the modification and appUca- 
tion of the Hodgldn41uxley equations to other excitai)lB tissues. In the 
hearty for exaihple, equations of the Hodgkin-Huxley^ type deacribe the 
electrical processes that generate the cardiac rhythm and^^ coordinate the 
heartbeat. Mathematicians are just beginning to use th%se equations as 
a basis for a theory of the abnormal rhythms of the b^eart, of which the 
most serious is ventricula* fibrilation. This theory has connections with 
recent Vifork on chaotic dynamical systems and transitions^ turbulence: it 
appears that fibrilation is directly analogous to turbulence. This work has 
enormous practical significance, since the principal cause of sudden death 
following heart attacks is ventricular fibrilation, 

Progress ,has also been' made in the modeling of the input to the 
neuron (whose output signal is the nerve impulse). The neuron integrates 
information received through a tri^ of dendrites in which the signal- 
ing mechanism is often described hy the^Jinear diffusion equation with 
leakage. Mathematical modeUng of the dendritic tree-has ^ad a substan- 
tial impact on experimental neurophysiology. One reason for this is that*^ 
dendrites too small to be penetrated with microelectrodes. Thus the 
neurophysiologist can only record voltage or inject current at the cell body 
and is forced to rely on the theory to indicate the sigi^ficance of these 
measurements with respect to activity in the dendritic tree. 

Some major successes of the theory are as follows: . ^ 

« 

1, Elucidation of the dramatic differences between effects of a synapse 
dose to the cell body and effects of a simUar synapse far oi^t'in the dendritic ' 
tree. ' ^ 

2, A possible explanation of the role of dendritic s^nnes in learning 
and memory* * 

3, Prediction of the existence of dendro-dendritic synapses based 
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, on a mathematical model of field potentials in tbe olfactoiy bulb. Such 
synapses, previously unheard of^ were subseqtiently found in electron 
microg^pbs. This work led to a fundamental new concept of* local in* 



without nerv^ impulses. Sucb processing is important in neural networks 
such as the retina. Indeed, an extensive mathematical theory of the retina 
has been developed, and this is another exciting: area'^of current research. 



2.12 ELECTRONIC COMPONENTS 



The design^ and fabrication of modem integrated circuits is a complex 
process. The number of devices that one can put on a chip (|epends on 
tbe size of the chip and bow small one can make its features, X^r tbe 
years, the largest increase in the number of device oh a j!hip has resulted 
from tjtte continuing reduction in feature size and with this a reduction in 
device size*^ 

(Consequently, process and design oigxncers bave bad continually to 
redesign the process steps and then recalculate the resulting device charac* , 
teristics to ensure good^^lctrical behavior* This bas' had to be computa* 
tional because tbe equations are mathematicallyln{^actab^e and a trial and 
error approach is prohibitively expensive and time'consumingr Moreover, 
experimental techniques tell us only ivhat happened, not wb/r Effective 
device design depend^'on determining both tbeivbat and the why by vary* 
ing the problem parameterslh the computational model. 

The mathematical models on- which- tbe^b^ry of semiconductor 
devices Vests are differential equations that describe the flow of current 
(holes and electrons) under the influence of electric fields. When feature 
sizes were large, the devices could be treated as though they consisted 
of plane surfaces and edge effects could be neglected. This allowed the 
development and successful use of one-dimensionai analytical models ob- 
tained by solving a stystem of three coupled^ nonlinear ordSnaiy differential 
equations. As device sizes shrunk, these models became more complicated 
and IMs accurate as edge efiects became more important. In very large- 
sqale ii^tegration where device dimensions have reached a few micrometers, 
tbese modeb are jio longer adequatj^, ,an^ the coupled, nonlinear^ partial 
differential equations must now be solved in two and^^ggaetm three 
dimensions. ThesB differential equations consist of a nonliiiW Pdeson 
equation that describes the potential of the electric field^aiiLd two non* 
linear transport equations that describe the motion of the holes and the 
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electrons. Moreover, the smatler sizes have made some physical effects 
important whereas they were previously neglected. This results in further 
complications to the basic' equatipns. 

Partial differential equations play an increasingly important role 
in simulating the fabricatbn processes. The transistors in a chip are 
formed by implaftting; certain dopantjoj>^nto selected are^ of the chip. 
Subsequent high- temperature proc^s^s, such as growing aa 6xide layer, 
will cause these atomic fmpurities to difibee. Their final distribution is an, 
important factor^in determining device characteristics. Insighta into thede 
fabrication processes are especially important for' increasing the yield of 
reliable devices, which is a critical factor in*their economic viability. 

Th^ design of the overaU^cuit to be placed on a chip leads^ large 
^stems of nonlinear differenflplvequations that need to be solved numeri- 
cally. Then the efficient layout of the circuit bn the chip. introduces com- 
binational ^Qdgraph theoretical problelms, 'which again pose formidable 
computatioSTjN^gblpms. * ,^ ' ^ 
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Computational mathematics and modeling is a relatively young science, 
one that is expanding rapidly. The success and importance of the field 
stems ftom the fact that its application provides the possibility of tack-^ 
ling sijpiuScantly more complex and di£9(^t problems than would othei^ 
wise be possif)te. Perhaps the greatest opportunity provided by a com^ 
putational approach is that it opens the wide realm qt stroi];gly non* 
lihgar phenomena to ^stematic, relatively accurate, and effid^t models 
ing, improving the chance that important phenomena can be isolated 
and an^ysed. Nonlinearities pervade nearly all aspects of applied mathe^ 
matica, and to a large extent these nonlinearities are responsible for the 
difficulties that are encountered in computational moAUng, Our purpose 
in this chapte^is to expbre tjie sqprce of so^e of ihese difficulties. In 
Chapter 4, we will discuss the direction (of some of the computational 
researchndeeded to resolve the problems. ^ 



3.1 DEGREES, OF FREEDOM * ^ ^ 

There are several reasons for the d^rees of freedom in a m<idel, and 
hence the size of computational problems^ to increase. One is that we 
attempt to increase the accuracy and complexity of our representations of 
the physica! conservation laws. Increased d^r«^s bf flreedom come either 
«f rem increasing the number of dependent or independent variables. Fbr ex* 
ample^ in chemica! models the nunjber S' dependent variables is increased 
by increasing the number of chemical species^considered. .An bbvious 
need forincreased independent variables comes ttom the need to r^pce^t 
phenomena in two and three spatial dimensions. However, even higher-^ 
dimensional p^lems arise when the independent variables are not the spa* 
tial cooi^dinates but are various state descriptors; such hig^r-^ditnensional 
proJ)lems are common in pli^sics and chemistry, Unfortunatelyftjie direct 
application of numerica! methods that work, well %£^e ojr two dimensions 
often are not. usably in three dimensions. Therefore/&Scmsi;ig the degrees 
of freedom may require significantly different algoritbqis:\r. 



Another way to in^j^ehse the degrees^offreQdom.in prohlems is to. 
ailow the possihility of multivalued adutions^^Certain physical Bystenls 
aUow the solutions^to hifurcatCj a phenomenon discussed in Section 3-5. 
Other systems have hysteresis or "bistory-dependent";properties. In these, 
problems the solutions depend hot only on the boundary conditions hut 
also on the path that the transient solution follows. Structural dynamics 
models for plastic materials in which the stress depends not only on the ^ 
strain 6ut on the straining rate as well exhiljit this behavior. 

Consider an example to quantify the magnitude of the ccsnputational 
requirements for multidimensional problems. Fir«t suppose we want to fol- 
low the evolution of a chemical model having N chemical species. Suppose 
also that the problem is stiif (see Section 3.2) so for each time step ia 
the evolving system we solve N nonlinear algebraic equations in N un?.^ 
Renowns. Even for models with ^ l^i^number of species this is a relatively 
straightforward task. However, if tVe now want to introduce transport 
phen<wena, such as Quid miKingf the problem has to include the spatial 
dep^endence of each chemical species. ' ^ ^ * 

For a one^dimensional case assume that yi§ use / spatial mesh.polnts 
and that w estimal^yrhat each species will do on the basis of its current 
localvalue and that of its immediate neighbors (a three-point jspatial stencil, 
involving all species). We now have _NI unknowns* and a nonlinear systetb. 
that has 3iV/ nonzero entries in each equation. It is typical to.have4pto50 . 
species and to require^OO m^K pointy to resolve the species c(^centrations, 
accurately. In this case over th^ quarters of a million words of memory 
are needed Just to store the approximating local linear system along with 
the solution. This alone is.Istger than the fast-access memory of most 
modem computers. Suppose /im;her that the same model is to be posed 
in two dimensions on a / X/ mesh apd in three dimensions on sXi Ixixl 
mesh. The solutions themselves require NP and NI^ words of storage,^ 
respectively. Worse, the totality of coefficients in the typical linear system 
will be 5N^I^ 6jid IN^I^. It is of practical interest to want solutions to 
three-dimensional combined Icinetics and transport .problems. However^ 
even for a modest system of 20 species and 50 mesh points per spatial 
dimension^ 2.5 million words are needed to store the soluftion alone. An^ 
additional 350 million are required to stpre the coefficients of the linear, 
system. (For comparison, we note that the largest computers. currenttjr 
available have .4 million words of fast-access rneiffory.) Mbreoverj^ iii the^ 
high'dimensional ca$e$> the linear system is not conveniently structure 
for efBcient solution. . / 

Cle^ly a major problem in solving high-dimensipnal ^steml^^of par-^ 
tia! differential equations is that after discretization the resulting.^tem. 
of approicimating linefu* equations can he much top large to he aoly^ . 
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effectivBly by direct means; Except, in simple cases the resulting equations _^ 
have to ^ soWed iteratively, and, especially for strongly nonlinear e(^a* 
tionsj the cpnvergence properties of ^e iterative process may be a major 
concern. ^ 

Operator splitting methods, such as alternating direction implicit 
(ADI) (see Section 4J), also provide a current effective approach for 
high*dimensional problems, T^e approach here is to alternate the sola* 
Uon of a set of lower-dimensional problems. Jhe alternated problem ap^ 
proximates the original problem vnth sufficient accuracy, but the set of 
lower*dimensional problems is much more, easily solved^ even in aggregate. 
Again convergence and accuracy properties must be established in all but 
the simplest of cases. 

Fbr some problems i^ may ]be more efficient to depart bom the^ con- 
ventional ideas of discretization on a iniesh network and consider instead 
mesh-free methods or (se^^Section&'4.8 and 4;d)^Mpnte Carlo methods, 
tfnknpwns^^ be change^l: instead of the^moimt of material at a g$ven 
location,;on^ ask for the amount of a^given wavelength in th!^.8p!uUon 
as a whoTe^ In Suid. mechanics, vort^ m^hpds could be a mor^ f^dent 
approach. Such methods reduce the size of the linear algebra problems, in 
comparison ,:with the m^es^orient^ methods^ jf thanew unknomis car^y^ 
'^ali t^^impOrtanWirfoma^^ ' - '^ ^ ' ' 

l>eiisipns abouV which mjeth<Kls are most efiective may depend; 
strojgly .on adyfuac^^and . changes in computer architecture. Since new-_ 
arcm^ectures, such as vector and parallel processors, a^ now evolving^, 
the numerical a^yst has to re^evaluat^his approaches periodically. Also,. 
^ considerations sucn as the rejatiye cost of memory, y^us central processor 
time c^ bea^ heaWy^Lon 



Inp^nnciple a model cai; contain impoi^^t length scales that range ^m 
the size of an ^l^l>o.the size of tthe unive^.,L) prac^ice^. however/ we 
linut the range of scales by f^pracimation (see> for example. Section ^ J). ^ 
^Nevertheless,, it^ is often, the case.^thM 'ma^ematical models of physicai 
processes are cKarfiu:Wri^d^tv the presence hi their solution 

of significant^ 4^^nt^ time and length scales. The jsol^tion t^ese ^ 
mpdels^wUI ImveTegions pr^trcngly. localized^ 1^^^ such as shock?, 
st^ep. fronts, jor^ther near disconthiti|tieS^ Therefore inoiportant topics of 
research in numerical analy^s, are the consideratiob of such (^umstances 



and the development of efficient theories and methodologies for Uieir com- 
putational solution. Indeed, we often find stxiations for which solutions 
aie not possi|;kle, or at least not practical^ without the application of spe^ 
cialized methods to deal with the multiple characteristic scales. 

Fluid mechanics and chemical kinetica are two areas that provide a 
rich source of examples for multiple anddisparate scales. Fluid- mechanical 
proc^ses are conunonly characterized by groups of par^eters that are 
indicative of the various scales ia a problem. Some examples are the Mach 
number (relates velocities to sound speed), the Reynolds number (relates 
inertial forces to viscous forces)^ the Prandtl number (relates viscous effects 
to thermal effects), and the D&mkobler number (relates chemical reaction 
rates to diffusion rates). Wb^ any of these numbers is very large (or 
small) it is like^ that'the solutions to the models wiU have regions of 
localised behavior. For the case of large Mach liUmber the possibility of 
shocks exists. Similarly for large Reynolds number we expect to ^counter 
boundary layers in the vicinity of solid walls. When the Prandtl number, 
is large we expect that thermal boundary layers will be much thinned 
than viscous boundary lag/fers. And when Damkohler numbers are large we 
expect narrow reaction fronts. 

In chemical kinetics we find large numbers of chemical reactions com- 
peting for the available chemical species at ividdy different rates. M a 
result some species are either consumed or produced' rapidly or slowly, 
white other species are being both produced and consumed simultaneously 
at high rateSt .with their net production rate beings reJativdy slow. Thij^ 
chemical behavior is responsible for the mai^ widely differing time scales 
in the mathematical models. The computational models of these processes 
are characterized as either multirat^'problems or stiff problems. 

It is worthwhile, to point out the distinction between, multirate. prob- 
lems and stiff problems. In both problems the system itself is- equally 
capable of rapid or slow changes. Multirate problems are those in which 
at least one' component of the solution is changing rapidly^ evm though 
others are changing slowly. Numerical methods for these problems n)^ust 
take time steps that are small enough to resolve the fast transients^ so they 
are controlled by accuracy not stability considerations. Stiff problems, on 
tlie other hand, are. those in which all components' of the solution are 
changing slowly compared with the fastest characWistic scales possible in , 
the model. In these cases ^plicit numerical methods are forced to take 
much smaller time steps than are needed to maintain accuracy in order 
to maintain stability. Oftep problems that l^gin as multirate problems 
become stiff problems as an equilibrium or steady-state condition is apr 
proached.^Stiff problems are usually solved efficiently by implicit methods. 
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For stiff or multirate problems, it is perli&ps useful to consider methods 
that treat the fast components differently ihm the slow components. 

Given that multiple^scale problems are of practical importance, we 
^^^must consider vthy are there computational difficulties, and what can be 
^aone to ameliorate the difficulties? the case of diBparate length scales 
one difficulty is that of representing accurately the highly localized be- 
havior. If. the solution is represented discretely on a mesh network, the 
mesh must be sufficiently fine to capture,the localized behavior accurately. 
The whole topic of adaptive meshing is critical^ important for these prob- 
lems (see Section 4.2). Here, instead of computing on a fixed prespecffied 
mesh, the mesh adjusts itself dynamically as the solution develops in order 
to maintain accuracy in the solution. ^ « 

For sit^tions in which the localized behavior is known to be ap- 
proximated well by veiy sharp fh>nts (e.g., shocks or flames), front track- 
ing methods can b&ve significant adviatages. Unlike the adaptive mesh 
approach ivhere the solution is resolved smoothly through the flront, the ^ 
front tracking method^ approsdmate the front by a discontinuify ivhose^ 
ma^tude, speed, and location are to be found. Then elsewhere in th^ 
^^fgion the conventional discrete representations are ade<}uate. * 

In problems like chemical kinetics, the disparate'.time scales cav 
the gov^ing differential e<}uations to be stiff Here, explicit soluf/ion 
methods are^well known to be extremely inefficient, and some for 
implicit method is needed. For systems of ordinary differential equ 
the problem haa been worked out, and bigb'-qualify computer softv 
available. However, %vhen stiffness is encountered in the context of i 
of partial differential e<}uations the remedies are much less developed, 
same techniques used for ordinary differential. e<}uation8, when applied 
directly the partial differential equation problems, of practical interest^ 
often yield problems that are simply too targe fbr current computers. 

Several approaches show promise. One is to develop pperato^spl^tting 
methods in which the stiff^arts of the problem are split off and solved as 
a series of^ smaller and h^ce more tractable problems. Another approach 
is ta attempt to remove the stiffiiess^by solving instead an approximate 
(yet sufficiently accurate) system of ecjuations. This tack benefits from an 
asymptotic analysis qt the equations.. Boih approaches have found recent 
successes in fluid mechanics and in combustion chemistry. 
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3.3 SINGULARITIES IN COEFFICIENTS, DATA, OR STATES 



Difficulties similar to tbose encountered in multiple-scale models are often 
found in problems having singularities In coeScients or states. That is, a 
singular or discontinuous coefficient can give rise to localised behavior in 
tbe solution, sucb as very steep fh)nt3. An example could be a material 
interface in a structural or beat-transfer problem^ say between a sted and a 
plastic part. Attbis material interface the solution (stress or temperature 
gradient) might change rapidly. In order to maintain accuracy in the 
computed ^lution, the numerical procedure would have to resolve this 
frontlike behavior. The situation is analogous to the occurrence of a shock 
or a flame front. However^ we usually Iqiovv where the material interfaces 
are, so they are perhaps easier to deal with than phenomena such as shocks. 

Phase transitions can also produce discontinuous coefficients. Take 
a model in which a melting fh)nt is traveling through a region. Usually 
the properties of the molten material are quite different than those of the 
soUd material. In fact, different governing equations may even berequired 
for the two regions (e.g., fluid motion may be .modeled in tbe liquid hut 
not the sohd). In any case, tbe solution is likely to exhibit jumps in its 
properties at the transition^ and tbe numerical method Avill need to locate 
and resolve it. This situation is more like a shocks in that the position of 
the phase transition front is not known a priorij and thus the numerical 
method must both locate and resolve it. 



3.4 BOUNDARY CONDITIONS 

The solution to a boundary, or initiaUbounjlary-value problem depends 
strongly on the boundary conditions. Thus it is important to understand 
the relationship of the boundary condi^ons to the differential equations 
and to tbeir discrete representation. Most important, the boundary con- 
ditions must be chosen so that tbe problem is well posed. For a large class 
of problems there is a satisfactQi:y theory of admissible boundary condi- 
tions, but for many problems^ those involving coupled fayperboh^eUiptic 
^stems or disparate time scales, for example/ only a rudimentary theory 
is available^ ^ ^ 

A possible error in prescribing boundary conditions for hyperbolic 
equations is to overspecify or underspecify the number of boundai^ condi- 
tions. Overspecification usoaUy causes nonsmooth solutions with mesh os- 
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cillations^near the boundary. Underspecificatioi^ does not ensure a unique' 
solution, and the numerical solution may tend to wander in steady-state 
calculations. In either case, the results axe not accurate and one should 
be skeptical of even the qualitathre behavior of the Bolution. It should 
be noted that the \vay in which boundary conditions are specked for the 
difference equations can change a well-posed continuous problem into an 
ill-posed (unstable) discrete problem. ./^ 

yTwo of the moat common methods used to incorporate boundary 

...^conditions into discrete equations are the extrapolation and uncentered 
differences methods. In the extrapolation method, the domain of the prob- 
lem is extended and the solution is extrapolated to fictitious points outside 
the integration r^on. The nonphjrsical solution at these points is defined 
so that the discrete equations are consistent with as many relationships 
as^ can be derived from the boundary conditions and differentia] equa- 
tions. The extrapolation formula can do this best by incorporating the 
discrete boundary conditions^ into the extrapolant. Additional relations 
can be generatedby differentiatibg the boundary conditions with respect to 
time, replacing all time derivatives by space derivatives using the original 

^ differential equation, and discretising the resulting eijuations. 

The uncentered differences ap|>roach is to extend the number of 
boundary conditions so that all components of the solution .axe defiiied 
at the boundary^ Again, these additional boundary conditions must be 
consistent with the original problem and as many relationships as csxk be 
derived from it- An uncentered difference approximation is then used to 
approximate the spatial derivatives at the mesh points nearest tKe bound- 
ary. V ■ _ " , 

Irregular domains can be imbedded in an underlying regular grid 
that is not aligned with the boundary, or.ari irregular grid can be con- 
structed that conforms to the boundary^ The discrete appr<^dmatipns to 
the equations away from the boundary are rpuch simpler on the regular 
imbedded grids, but the boundary coiuiitions are difficult to approximate. 
Boundary-fitted grids can be generated algebraically in the original ph3rai- 
caj domain^ or the domain (and hence tlie grid) cm be mapped onto a 
regular tgrid in a simpler domain jnd the equationa solved^tbere. The 
algebraiC'gfidTgener^tion mett^odshave the advantage that the equations 
and boundary cqndiUons are unchanged^.but the dtffereijitial bper^torsare 
more diiSciuilt to .approximate on the nonuniform grid. When using the 
mapping method, ttie diETerential operatio^is are eerily approximate but 
the transforinatioit can.|;reatly compUcate^the equation and sometimes 
obscure important- prop^ies such as'tKe conservation la^ expressed by 
the equatiohs.^^: ^ \:; ■ . 



3.4.1 Boundary Conditions at Infinity 



Many physical problems require the solution of paitial dtiferential equa- 
tion^ on some infinitely large domain 0. For computational reasons this 
domain is often replaced by a finite domain Oi- Then the difficult problem 

^of ^ecifying boundary conditions at its finite artificial boundary B arisea. 

^ti^^especially important that these artificial boundary conditions do not 
Introduce spurious phenomena* Consider, for example, a nonviscous iluid 
that at Subsonic speed leaves Oi through the boundary B. There, is one 
characteristic direction that points back into the region 0^ and there* 
fore one boundary condition has ^to be specified on B. But, in general, 
no detailed knowledge of the solution on B.is known and therefore other 
principles have to be applied. For example^ if one has solved the problem 
by diSerence approximation then one predicts the solution on B from in* 
side Wit for ail the dependent variables, ^hus this procedure amounts to 
overspecification of the solution on B. Another prinjnple has been proposed, 
namely^ to specify the boundary conditions on B so that no reflection of 
high frequency take* place. Howe^r, numerical experiments have shown 
that suoh approaches ido not always work. 



3.5 BIFURCATIONS AND CHAOS . ^ 

3.5.1 Buckling and Collapse Behavior, Bifurcations 

In general^^the equilibrium equations of a mechanical structure involve a 
finite number of pSfameters, that is^ they ime the generic form F{xyp) ^ 
0, where x varies in some state space X and p 6 it^ represents a parameter 
vector. Thas, in general^ the solution set {(a;,p)"€ X xR^i J^fej?)f= 0} is 
a manifold JCxiif^^ and one topic of interest is the location and chai:acter. 
of the singiilar points of'lhe'solutiohlet.. ^ ^ " > ^ 

Without going into detail^ suppose that we follow some curve on 
the solution manifold defined by some combination of parameter values^ ^ 
with one degree of .fr€?dom represented by a single variable X. Then - 
we may encounter certain critical points on the solution curve where 
the mechanical stnictu^ may suffer a loss of stabfiity. Such a loss of 
stability actually corresponds to a dynamic phenomenon whereby the 
stnictureundergoesasuddenchwgeof defonpation. The dynamics of this 



phenomenon are not described by the equatiotis of equilibrium F{x, p) ^ 0, 
but it is possible to deduce from the shape of the (static) ^lution manifold 
at that p(»nt in4iat type of sudden changes may b^ expected. 

We use a few figures to illustrate the situation. In the upper left-hand 
part of Fig. 3,1 the point denoted by (1) is a so-called limit point or turning 
point and att increase of the load-intensi^ X beyond the critical value Xc 
^results in a jump from,(l] to (presumably) (2), This type of behavior is 
called a snap-through or collapse. In the case of the upper right-hand part 
of Fig, 3,1 the instabili^ phenomenon is related to the bifurcation of the 
solutions at (1), This behavior is called buckling. This pe^ of the figure 
is a classical example of stable buckling where a distinct change in the 
character of the load deformation is encountered when the toad-intensi^ 
X passes the buckling load X^ but where no failure of the structure occurs. 
On the other hand, the lower part of Fig, 3.1 shows an unstable buckling 
point where we again observe a dynamic departure from the Ijifurcation 
point to some other state [presumably the equilibrium state (2)], The 
geometrical shape of the bifurcation branch E is the determining factor in 
the question of whether the bifurcation point is stable or not. 

These examples already indicate'that for a deeper understanding of 
the behavior of a mechanical system it is necessary to analyze the sEape of 
the full solution manifold. Of course, the choice of the parameters entering 
into the definition of the equation is of critical importance here. In essence, 
catastrophe theory provides sonie information about the selection of par- 
ticular minimal numbers of parameters, but in practice the parameters 
are ' simply chosen to correspond to the certain natural features of the 
structural problems. 

In view of these observations the aim is to develop procedures for 
a computational analysis of the Torm of the equilibrium surface. Some 
methods for this purpose ar^ mentioned in Section 4,4 on continuation 
methods, but the entire problem is still a wide open research problem, In"^ 
particular, it has to be noted that we can compute only points that are 
approximately on the solution mamfold of a discretization of the original^ 
problem. Thus, the questions arise whether phenomena, such as limit 
points, or stable or unstable bifurcation points, encountered on the solution 
s^t of the discretized problem actually correspond to similar phenomena 
for the original problem, and, if so, what errors have been encountered. 
These questions are as yet l^^eljLun^swered and represent considerable 
research challenge. In this connection, it might be mentioned that for 
nonlinear problems of this type the solution manifold of the discretized 
equations often has a different number of connoted components than 
that of the original equations. The componenj^ that do not approximaU 
the exact solution manifold have been called spujious^ or numerically 
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FIGURE 34 Instability behavior of shells under different loadings F, On 
- the right the load intensity X is plotted versus a chwactedstic deflection, 
(fVom W- Rheinboidt and E. Riks, ''A sumy.;o{LsoUd techniques 
for nonlinear finite element equations," Si(xit of the Art Surveys of Finite 
Element Technotogy,KC^Bpr 16; speci^J AMD publ., American Society of. .11 
Mechanical Engineers, New York (1983).] ' 



irrelevant solutions. They are being observed more and more often in 
various applications, but their study has only recently begun. 



3.5.2 Chaos in Deterministic Systems 

It has been remarked that one of the most engaging problems in nonlinear 
dynamics is that of understanding how simple deterministic equations can 
yield apparently random solutions. As is now widely recognized, systems 
with this behavior appear in many of the branches of scientific endeavor. 
This near universG^ity has come to the attention of physical scientists and 
applied mathematicians principally through the use of computers in the 
study of properties of dynamical systems. 

The onset of chaos in' deterministic systems (the s^hastic instability) 
signals an unusual sensitivity to initial conditions^ i.e., the trajectories in 
phase space diverge widely as time goes on even though the initial con- 
ditions are arbitrarily dose to one another. The behavior is such as one 
would expect in a space with everywhere negative Gaus^an curvature. On 
the other hand the trajectories may in fact tend to a single orbit nearly 
filling a subspace of the phase space (a strange attractor) in an ergodic 
manner. In even the simplest systems of this type it is possible for several 
such attractors to coexist side by side, with the initial conditions deter- 
mining which one is reaittied a^mptoticaJly in time. These and other 
numerous unusual properties of dynamical systems displaying chaotic be^ 
havior have led to some oAen problems in computational mathematics. 

Examples of such problems include the following: When on vaiying a 
param^r the stochastic Stability sets in, the continued use of an algo- 
rithm describing the evolution of the system prior to the onset of chaotic 
, behavior may no longer beappcopriate. In this event it may be more piac'* 
tical to'tal^advantag^ of the stochastic nature of the ^atem and use more 
or less conventional statistical methods. In order to eSect such a change in 
computational methods it is necessary to detect the change from orderly 
behaivior of the ^tem to chaotic behavior. ^Theoretically ^ the sensitivity 
,or the system to initial c<»iditions accompanies the occurrence of positive 
Lyapounov exponents," i.e., numerical indices of the a^yniptotic rate of 
divergence of initially nearby system trajectories. Since the Lyapounqy 
exponents are defined as time asymptotic quantities, straightforward com- 
putation of these quantities requires a.considera)>le computational effort, 
including a calculation of the evolution of two initially nearby tr^ectories, 
or the simultaneous integration of the associated variational equations. 

An important toot for studying the properties of dynamical systems 

57 , 

ERIC : , ' ■ • > . 



in the chaotic regime is the first return mapt sometimes referred to as 
the Pqfncar^ map. Such a map displays in a graphical manner various im- 
portant aspects of the attractor. Unfortunately, the BtmpUcity of generate 
ing such a map is restricted to systems no larger than three dimensioDal, 
because of the self-evident difficulty of making multidimensional graphic 
displays. As the study of higher-dimensional dynamical ^stems advances 
there wiUJ>e an. ever more urgent need for-a higher-dimensional equivalent 
of the first return map. 

As motioned earlier, in the chaotic regime it may be more appropriate 
to describe the properties of the ^tem trajectoiy in statistical terms than 
m terms of a trajectoiy evolution. In order to obtain such a statistical 
description it is necessary to have appropriate information about the in^^ 
variant measure associated with the given dynamical system. With few 
exceptions such measures cannot be derived a priori but must be obtained 
from detailed calculations of the trajectories. A finite machine computa^ 
tion of the ^stem trajectory; will inevitably introduce some errors that 
may be all the more serious in the chaotic regime because of the previously 
mentioned exponential divergence of nearby trajectories. There is then the 
qfuestion of how accurate the computation of a chaotic tre^ectoiy must 
be in order to yield enough information for constructing the appropriate 
^variant measure. ^ 

While the above examples have been drawn from dynamical ^s^ 
terns representsble by ordinary differential equations there are other sys- 
temsj e.g., those described by partial differential equations and by integral 
^ equations, whose study is likely to be replete with similar computational 
difficulties. 



3.5.3 Symmetiy Breaking 

Bifurcatioiis commonly arise in connection with a loss, or breaking of 
symmetry. In such cases the extra structure of the problem symmetry 
may simplify the analysis. We give an example from fiuid mechanics. A 
classic experiment concerns flow between two rotating vertical flinders. 
For small angular velocities of the (^Dnders the flow is laminar and is called . 
Couette flow. As {he anguUr velocity is iniueased the vertical translation^ 
symmetiy is broken and axisymmetric Taylor cells appear. These cells 
resemble a riow of adj^nt smoke rings, each one rotating in a sense 
opposite to its immediate neigh|)or. A fbrther increase in velocity, breaks 
the axial symmetry ^ and the Taylor cells become wavy with a time periodic 
shape. Eventually the Bow becomes turbulent. Specific experiments on 
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these systems are a classic topic in fluid meclianics. However as often 
' happens in science, one cannot see important phenomena in the absence o 
a cflrefully thought out theory. Recently a number of precise estperiments 
hav^ been conducted concerning the role of strange attractors. Routes to 
turbulence and, in particular^ wavy Taylor cells were carefully observed. 
The results confinned some of the theoretical predictions but not all of 
them. One theoretical proposal was that generic limiting sets, known as 
strange attractors^ would dominate turbulent flow.pattems. This proposal 
was supported by the fact that strange attractors were observed in low* 
mode approximations to fluid flows. More refined calculations with more 
modes included indicate an absence of these strange attractors, but the 
phenomenon is still indicative of the highly complex solution mamfolds 
that can arise in nonlinear problems. 



3,6 ILL-POSED INVERSE PROBLEMS • ' 

r * 
The notion of a well-posed^problem is due to Hadamard: a solution must 
exist, be unique, and d^^d ct^ntinuously on the data* The tma '^data" 
can have a variety of meanings; in a differential equation it could include 
any or all of the following: boundary values^ initial values^ farcing term^ 
and even the coefficients in the differential equation.- Since data cannot be 
known or measured with arbitrary precision^ it was felt for along time that 
real physical phenomena had to^be modeled by mathematically well*posed 
problems. This attitude has changed cotisiderably in recent years, and it 
is now recognized that many applied problems are iUr^sed^ particularly 
when they require numerical answers in the presence of contamination of 
the data, 

ni-posed problems often arise in the inversion of well-posed problepas. 
Consider, for instance,va well-posed problem that smooths the data or 
attenuates its high f^quencies. The inverse problem, in which the role of 
. data and solution are interchange, will then be ill^posed, A simple but 
important example is the F^%dholm integral equation of the first kind ..^ 

^ ' ft(x,y)u{v)^= v{x)^ 0 < X < 1 

6r, in operator form / * 

Assuming the kernel to be continuous on- the closed unit sqioare, the 
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Riemann-Lebesgue lemma gives 
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, liiri f k{x,y)^nnydy = Q 

JO 4 



SO ttiat K attenuates High frequencies and thus transforms widely different 
functions u into approximately the same v- is regarded as the data, the 
solution u of the integral equation, therefore, cannot depend continuously 
on the data. Indeed, it is not necessarily uniquely solvable, eitti&-take, 
for example, 

JTp unaware of these diflicullies, one attempted to solve the integral equa- 
*on by discretisation, one would find that the corre^ndin*g matrix prob* 
lem is singular or ill-conditioned (and the finer the discretization, the more 
ill-conditioned the matrix problem). 

« More general inversion problems can sJlso be reduced to an equation 
of the form Ku = v, where K is a eontmuous transformation, v is the 
data, snd u is the solution, being sought, m may happei^tha^ u itself is 
not the quantity of principal interest but rat^r some funcfionals of ti such 
as some of its momenta or its values at a fevr^gecified points. ' 

Thus, there are three pieces of information that are central to the 
numerical resolution of an inverse problem: (1) the model M, i^preseating 
the equation involving a mapping K between appropriate ^aces; (2) the 
observation operator 0, representing the measurements that can be made; 
for instance, we might have p{t;) = {t;(xt);i = 1,2,., .,n}; (3) the tnteU 
Hgence opetator J, which specifies the information^^we wish to extract from 
the solution. For instance, we might have J(ti)^ {/x*ti(a;)dicj fc=J[^,l,2}. 

Most of the existing approaches focus almost exclusively on the model 
M, taking account of a priori information about the solution such as 
smoothness, positivity, and bandwidth. The Tychonov regularization 
method is*of this kind. By restricting K to a suitable subspace it may^ 
be possible for the restriction of K to be one to one with a continuous in- 
verse. It might be more appropriate to study the tripletJVf|0,J for different 
choices of 0 and J. Ope gqal of such a n'ew approach would be to classify 
and quaittify ill-posedneas by developing comparison principles and prder 
structure, ' . ^ . 

Ql-ppsed problems with their attendant numerical difficulties abound 
in practice. For instance, in scattering from.obstacle^ one may wish to 
determine the^shape of an object (or its surface impedance, if its shape is 
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known) ftom far-field measuremenU. Coefficient identiflcation problems 
arise in many contexta, one such problem being the determination of the 
sound speed in subsuiface media from measurements of the field at or 
near the surface. Another class of ill-posed problems is associated with 
image reconstruction (tomography the medical applications of which are 
now widely used (CAT scanners, for example) (see Section 2.12). More 
recently, similar ideas have been applied to the nondestructive testing of 
mechanical structures to detect cracks in fuel rods, weaknesses in rotor 
blades of jet engines, and faults in screws in ships. Optimum filtering 
and inverse problems in Fourier optics (restoring data designed by a band- 
Umited filter, for example) are other areaa of current reseafth. 



3.7 EFFECTIVE MEDIA 



The study of "bulk** or ^effective** parameters for coi^iposite media is of 
fimdamental importance. Depending on the particular appUcatjon area it 
may be relevant to consider a periodic (or an almost periodic) or a random 
formulation. ' 



3.7.1 Homogenization, the Deterministic Approach 

Many probl^s of physical interest involve several length scales. As an 
important example we mention the study of composite mat^als in stnic^ 
tural mech^cs. Owing^to a particular manufacturing process, a distinct 
structure, e.g., periodically or almost periodicity, is often present. 

Homogenization i^ an approach for deriving ^ the . macroscopic 
properties of the material ftom the known microjscopic ones. A variation^ 
of this is to replftc^ a complicated geometrical configuration by asimpler 
one, e.g./to replace a framework by a plate or to smooth out a rough 
surface. " ' ^ * 

Homogenizatipn may be applied to Unear and nonlinear problems and 
, :^can provide, qualitative information about phyfflcal macrolj^vs as well as 
:^^^od approximations to the various parameters appearing in theyse laws. 
l^. '^he most important mathematical tool for thadeterministic approach 
is some form of asyiJipCotic expansion. The'tfieotetical results concern the 
limiting behavior when certain of the relevant scales become very small. 
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One ia likely to obf^n difierent homogenized formulatioDS from different 
limiting relations between the length scales^ 

The numerical solution of a problem ictroduces^i' new Jength scale: 
the mesh me m a fiiy^element method, the step size in a fimte-^difference 
method, or the wavelength m a spectral method. Different limiting rela-^ 
tions between this new length scale and the original.pbysical length scales 
Usually lead to different results. In t|iis sense the numerical treatment has 
to be considered simultaneously with thehomog^nhation process. One has 
to design algorithms that will adaptively select the correct homogenised 
formulation and discretize it appropriately. 

Whereas much theoretical woric has been done within the last ten 
years ana^sing the eEfects of multiple scales in continuous media^ the 
study of numerical discretizations of such pr<)blems is still very muclf in 
its inf^cy. * 



<H,i Analysis of Random Media 

Assessing' the efTeot of random fluctuations in the coefficients of a pa^^ 
tial differentia] equation is a basic mathematical problem that arises con- 
stantly in science^ applied science^ and engineering. Findlug effective con- 
ductivities of composite conducing materials such as^il or a metal alloy^ 
finding effective fluid equations for flow in porous media (Darcy^s law)^ 
determining the rate of sedimentation of particles in a fluids and many 
other pr<)blems are pr<)blems that may require analysis of random media. 
In addition to the quaiititative aspects of the problem^ many interesting 
qualitative questions can be posed as well. For example, what is the nature 
ot the spectrum of the Schrodinger e<)uation with a random potentialf^It 
has been shown that in one dimension^ for a large class of raiulom poten- 
tials, the spectrum is always discrete. For randomness of large magnitude^ 
it was recently shown tKat there is no difhsiour without restriction on the 
dimension of the space. 

From the appUcatlqns^ viewpoint^ one can frequently model ade- 
quately a random medium^ such as a suspension^ for example by a con- 
tinuum with suitable constitutive properties. Ttie continuum equations 
may be linear or nonlinear, and the constitutivg^^aws may be known 
only qualitatively itom experimental data. In such a context, theoreti'^ 
cal investigations are useful ^d. necess^^in order to understtod how 
the phenomenological continuum equations arise from the known micn> 
scopic structure. One can then fiiid mathematical characterizations for 
the relevant constitutive laws that can lead td interesting conclusions^ 



Examples of this arise fVequent^ as, for ^mstancei the detennination of| 
bounds for the eGfective conductivity of composites* 

The eCTectivg medi^ ideas of Mflxwell m the last ceotiuy domiaated 
theoretical calculations for & long time ^ In the laat two decades^ press- 
ing teclmological needs have caused a m^or expansion in materids re* 
search. In attempting, to achieve a mathematics understanding of the 
methodological basis for efifective media calculatioi^s^that abound^ one finds 
a lack of theoretical foundations. * 

Mathematical methods in random^ media are drawn from analysis, 
probability theory, asymptotic methods^ and di&rential equations. The 
goal is*to develop tools that bridge the gap between microscopic and 
macroscopic descriptions, give qualitative information about constitutive 
laws, and determine when residual stochastic effeeta remain and how they' 
can be characterized. In addition, it js of interest to find important specific 
problems on which more derailed analysis, .including numerical analysis, 
can be earned out so that these problems can serve as benchmarks in the 
theoretical, development. 



3.8 VALIDATION, ERROR ASSESSMENTS, AND 
SENSITIVITY ANALYSIS. 

The results look good, but are th^ real^ right? The question is often 
not answered satisfactorily, ^deed, often too little attenUon is paid to 
the difficult topic of model and code validatioiL Once the results look 
plausible^ we are often either unable or unprepared to take the valida* 
tion process (brthef. The validity of a model dep^ds both on having a 
proper pt^sical model (Do the governing equations adaquately repres^t 
the physics?) and on having an accurate computational representation of 
the govemipg equations. The mathematician h^ a responsibUity in both 
areas^ He should help determine the Vell'^posedness*' of the. models and^ 
ftom a mathematical point of view help the physical sd^tist determine 
the appropriateness of the model. Given that the equations are proper, the 
computational mathematician must be sure that the numerical procedures 
used accurately approximate the solution's|^ of the govjerning equations. 

FVequent^ modeling is done in coi^unction with expeBments, and 
those results are used to validate the model. Coc^paii^n with experiments 
'simultaneous^ testa both the validity. of the governing equatiojis in the . 
modet and of th^ numerical solution. In this case the source of any 
disf^pancy is not easily isolated. A reasonable gq^^should.be to validate 
the numerical procedures and the physical models separately so that, a 
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model can be developed more confidently. Theideaign of reliable numerical 
error estimates for^e computational methods is an impojrtant step. 

3.8.1 A Posteriori Error Estimates 

An important current research topic is the development and application 
of a posiertori error estimates. Such procedures are valuable for models in 
which differential equations are approximated by their discrete analogs on 
a mesh network. Wien (^posteriori estimates of the error associated with 
the discretization of the continuous model are ^ssible— and they often 
are— then potentially it is possible io control the errors, Adaptively moving 
the ^esh network to control or miminiiae the error is one application of 
this approach. With a posteriori error estimates, it should be possible, 
in principle, to give a strict bound for the error on completion of f£e. 
computation. 

The selection of the speciBc accurac3r requirements depends strongly 
on the goal of the computation. Often it is desired to obtain detailed 
information about fte solution itself^ in other cases, the main focus is the 
value of a specified functional of the solution, as, for examf^, a stress in- 
tensity factor in fracture mechanics or a drag coefiScient in Suid dynamics. 
Other goals may be the determination of certain critical data, such as 
collapse points or buckling points and their associated loads in structural 
mechanics. Sometimes special tecluiiques; $*g.^ variational techniques^ can 
-^]^eld highly accurate estimates of psfticular fiimctionals. 

In connection with most of these goals we are interested in quantitar 
tive results that have a desired accura<^. Fbr this the error has to be 
defined^ in that a family of exact rosults has to be specified with which 
the computed data are to be compared, a norm has to be prescribed in 
terms of which the error is to be measured, and some pr<ioed}}x^ has to be 
established for eatjpiating the error. 

Such error estimation capabilities are certainly importjuit in manv 
applications, provided they can be guaranteedlto be reliable. For example 
for the many types of certificatjon computations required in £He design o( 
complex structures or nuclear plants^ the availability of reliable estimates 
of the accuracy of the computed data is obviously important. In other 
cases, such estimates may reduce the total design effort and avoid un- 
necessary dverdeagn. At the same time, the availability of effective er^ 
ror estimates introduces the possibility of applying adaptive techniques to 
structure the computatioiji to achieve a desn^ed error tolerance at minimal 
cost or to provide the l^t possible solution within an allowable cpst range- 
Many of th^ theoretical studies of solution algorithms for classes of 



differential equations provide for some a priori error estimates, usually 
of an a^rmptoti(Lt^rpe. For the practical error assessment these a priori 
bounds are rarely compu^ble or Very accurate* Thus, one is led to the 
necessity of developing a posteriori error ^tin)atesl;hat utilize information 
gained duriBg tlM&coinputationutself.- ^ 

Fpr .th^ finite^ement solution of certain daases of elliptic boundary* 
yatue '^i^temSf-some computa'ble and reliable a posteriori error estimates 
have bee^^ey!dop^3^^^yzed in recent years. Most of these a po$* 
feHon estimates are ^asedo^^^ local an^ysts. For a given mesh A consist- 
ing of elements Ai,i..^An and Ayith a corresponding approximate solution 
on that mesh, an error indicator is associated with the jth element Aj. 
These have to be computable in terms of information about the prob- 
lem and the appro^dmate solution on Aj and,^t most^ on the Immediate 
neighbors of-that element. On the basi^ of the indicators riij^^^,r}n an 
error estinti^AfA) is then confitnicted. Of course, the Vj'as well as e(A) 
depend on. thejchosen norm. The theoretically important question is then 
the relationship between e(A) and the norm of the actual error e(A), The 
eEfectivity dt the estimation may be judged by the effectivity indes & = 
e(A)/t}e(A){|r In practice, it is usually more important fop^ to be^close 
to 1 tli^ that 0 <. 1.. Moreover^ it is essential that 0 converges to 1 when 
||e(A)|| tends, to 2er0r so that for an accuracy of^ BBy^ 5 talO percent the 
vaiue of 1^ 1| is j£xpected.tp-belessi;han| say^ 0,1 br 0.2. 

The {^ecfbtirep^'nt of designing error estimators with these properties 
for realistic blass^'of proldems and various different norms certainly repre- 
sents a demlinding research ta^k. The results' available so far suggest that 
such error estimfi^TB canl^ developed at least for linear proUems. For 
nonlinear probliemrthe situation is in an embryonic stage although some 
results available fcj model problems in one space dimension Indicate that 
estimators for the ^for'along contmuation paths and for thelocation of 
critical points aret^^^^'tstionally feasible, ^here is certain^ considerable 
need for concentrat^d^r^search in the.general area. 

3.8-2 Other Validation^ieasurer . * 

Checking for mesh dependencies ^nd^conv^gence rates is another- w^ to 
h^Ip validate a modeL If a method is hippj;^5ed to.be/say, ^second order ^" 
then as the mesh intervals are halved,, thfe^^roi^ by a 

factor of 4. i^'this does not hagl^ a$^»rdmg to one's error estimate^ then 
the method does not bave the desired order yet and may be in error. Also 
one should continue to refine the meish until the resulta^are independent of 
refinement to wi£5tl\ the desired accuracy. Solutions can als^ be checked. 



.for aensitivily to mesh orientation. The solution^ should not be aligning 
itself with some special property of the mesh network. 

Numerical ^nethods can also be validated by compahng.them ^th 
analytic aolutions. Often this is possible by choosing limiting cases where^ 
an asymptotic analysis applies. In some cases complex constitutive laws 
can be replaced witti constants so an exact solution is possible^ Higher-., 
dimensional models can sometimes be validated by comparison vdth pr^^ 
viously confinned lower-dimensional results. Bounds on solutions are 
another way to help verify a model. Often conservation laws or minimum 
and maximum principles apply, or the solution is known to approach a 
previously known equilibrium or steady-state condition. These properties 
can be used as a basis for validation. 

To the extent possible one should always write computer codes in 
a modular fashion such that each part of the^model can be validated ' 
separately. Even though this is just good structured programming prac- 
ticCf and not necessarily related to mathematics, it is an important aspect 
of model validation. 



3.8.3 Sensitivity Analysis ^ . ^ " . 

In addition to the independent and dependent variables, most models also 
depend on certain physical parameters. Unfortunately, those parameters 
are often not known accurately. Sensitivity analysis is a systematic 
means, to quantify the relationship of model paramete.rs and jaod^ results. 
Doing the sensitivity analysis requires solution of an ^additional system 
of differential equations. These equations are a formal Statement of 
the relationship between the dependent variables of 1^ system and the 
par-ameters. The results arfr^ven in terms of a itiatrix of partial deriva- 
tives of,the dependent variables with respect to the system parameters. 
One ^y to think of the analysis ia as a method to i^nrnde theoretical 
"errcff bars" for the model. " \\ ^ ^ 

New methods are being developed to solve these senStivity^equations 
quitelefiSciently/ The methods rely on,the (Jbservation that the ^sitivity^ 
equations are linear, regardless of the nonlinearities in the original nfodel/ 
The method sdves differential equations for the Green's function ^f^J^e 
sensitivity equations, and then the sensitivity coefficient matrices ar^ 
puted by quadrature for the various imhomogenous terms corresponding^ 
to the parameters. These methods have been successfully apphed to proK 
leins in chemical kinetics. 
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4. NUMERICAL METHODS 



4.1 DISCRETIZATION METHODS 

Many important pbys^ problems are modeled by boundary or initial;^ 
boundaiy-value problems. In this modeling, the physical state under 
consideration is characterized by a function Uy which is the unique solution 
to the boundary or initial-boundary'value problejn^^Thus a major part of 
the analysis of the physical problem is the detennl^tion of u or some of 
its derivatives or some functional.of it:^ 

In nearly all important problems the deternunatiDn of.u is an inSmt^ ,^ 
dimensional probleip in the sense that ^^d<^^ -not lie^iii an j^plicitly 
known finite-dimensional sjm or^ altern^ivety/cajo^nqt^^ in 
terms of a set of ekplici^l^cnawh function? by^iqeans^or^iinite set of 
parameters. Thus cannot be "completely** regimented on'a computer, 
and it is necessary toVesort tosomej^e of ap^roodmation or disctietiaation 
method. InessenfMyaUdiscretiaaWbnmetho^s^neattem^^ 
a function, Uopprox^X^^^.h^o'i theprys^Jiand is *c1q8C? to^t^iaAjon the other 
is ciiaracterized1by.a*finite.set of parameters, ini$o cah^be^culated^ 

The ^m of finding Uappfox.^that is "ctoS|S^ .^p ulis-Siaae precise by- ' 
requiring tftiat \\u - u^ppfox tl < t {or that \\u - Uappwix 11 < I'lt^tlDf where 

is a pryaically relevant norm and r is a physically relevant tolerance. 
The goal is^hus to find t*approx satisfying tWs grrol^cn^oii^,^^ _ ; 
least expenditure ,Gf^6mpjMtaU^^ * . ' 

procedure is infli^i^,^ ^uAibet^.^.considersd;t^s^^ unpd^ptant^ 
of which we novs^lii^*;^^ "7^/^^^.^ , : ^^^-^^^i^- ^^^i 



L fhe^goals of|hfeaii4ys|/<?M^e:p'fi^^ . i 

2. t^6wa.iaati^attcal p^j>erties ofjtbe physicsitiyrbble^^^ the 
algorithm:- ^ \ " \- ^ v,; X 4v'r r-Ci; A">.^ . ~'?P"^- . 

3. Hardw^e coi^side^ibn&i e.g.f the av^abiltty qf-pmlTel process. 

4. Uomputer^science considerations, e.g.7 data-management rec^uire- 
ments. * 



5* Restrictions on computation time and expense. - . 

► ^ - » , - * \ 
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^ As indicated above, the exact solution is approximated by a function 
t^apprax t which is expressible in terms of a set of explicitly known (basis) 
functions by means of a finite set of parameters, v^hich are determined 
Tn the computation. ^One brgad classification of methods is in terms 
of the nfLture oT the l>asis functioQS, namely (a) those involving global 
basis functions having global support and (b) those involving local basis 
functions having sms^ support. Another classification is based on the 
extent to which the method is adaptive. Adfiptlve methods refine and 
modify themselves on the basis pf partially complete computations. 

Of course, thfese classifications are not precise, and there are inethods 
possessing the different properties in various degrees. We now turn to a' 
more detailed discussion of the most commonly used discretization proce- 
duTes.* 



4.1.1 Finite Differences^ l^^^-^"^-^;- ^ rXjT 

_ . ^ . One of the most^freqyently used discretization nj.€thods is the method 
- ^ , of finite differences. Ttte .central ideajs to rej^lace each partial deriva- 
. tivei occurring in the differential equation in the luiderlying boundary or 
. Initial- bound^;^ue problem by an approximating difference quotient. 
. ^ For example, tha firstr order .derivative du{x,y)(dx may be replaced by 
. " the forwatd difference [^6'!^ + or the backward difTerence 

[u{x^y) - u{x - hty)]Jh, whereas the Laplacian An = Uxx + ^^3,3, may be 
» replaced by the five-point difference operator 

. ^ u(x- +h,y) + u{x, y + h) + u{x- h, y) + u{x, y - ft) - 4u{x, y) 
^" ^ _ ^ 

^This replacement of all derivatives^un the differential equation by ap- 
propriate difference quotients Je^s.^ i' ^tet? of.^ equatwns, called 
difference equations^ i^or the^it^mb.^'ju^, which axe to ^pl*c3anlate Jthe 
values of the exact solutipn^ tt(ar,jf),at the fimte^difierence'jiie^K points 
{jhfk}ij,j^i^fc = 0t^ij±2j^-.j where is a.sinaUpo^jve ntimb^ b^led the 
, m^sfi parameter., N^jmiform meshes may be ponadered as .wejS, tiiit at 
, cpos^crable compUs^fel^jn the form of the difference, quotients in.the 
. - resulting diffefence^"^'*'^'^- - " " '^^ '''^ ^' 



There are a numE^ of impoftanj; questions that arise iii c^^^ectipn 
with the study of difference metiiods, namely questions of . ^ 
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1. Accuracy of the method; 

2. Stability of the method; 

3. Efficient solvability of the finite-difference equations; 

4. Robustness of the method with respect to the input data (e.g., 
coefficients, forcing functions, meshes); and 

5. For nonlinear problems, the separation of legitimate approximate 
solutions, which correspond to an exact solution^ f^m spurious ap- 
proximate solutions, which do not correspond to any exact solution. 



Finite-difTerenc^metbqds are an example of the local basis class of 
approximation methods. They can be adaptive* Adaptivity is usually 
introduced via adaptive mesh selection: the mesh chosen at any stage in 
the computation is based on the previous computations. 

Much important work on difference methods remains to be done. We 
mention in particular the construction of effective difference methods for 
nonlinear problems, in particular foiuthoae problems whose solutions have 
^ocks. 



4.1.2 Variational Methods of Discretization 

We will now discuss a class of discretization methods based on variational 
Tormulations of the physical problem under consideration. As noted above, 
the problem of calculating the exact solutign u is iniinite dimensional in 
the sense that u is o^y known a priori Vf^e in an infinite dimensional 
space, S8y H. For a linear problem, a variational method of discretization 
consists of (a) a finite^dimensional space S c H, called the trial space, 
in which the approocimate solution is sought^(b) a finite dimensional test 
space V', and (c) a biUnear form £(u,^) defined onHxV'. The approximate 
solution is tb^n determined by requiring that 

Bin approx t ^) " B{Ui v)j for all V'€ V 

where for v€V, fl(u,v) isr'computable from the data of tbeproblem without 
knowing t£. Sin^S^ S and are finite dimensional, u^pptox can be calculated 
by means of the solution of a system, of linear equations if the o^i^nal 
equation is lifieaf! Usually^ 7or nonlinear problems^ B is nonlin^r^in u, 
and the resulting system of equation^ becomes nonlinear. 
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The approximate solution depends on the selectiofi of SiK^&nd B. 
For any problem there exists a wide variety of variational methods of 
discretization, i.e., a wide variety of choices for S, V, and S. The rational 
.selection of 5, V, aiid B is a central problem."' 

Variational methods can be of either the !oca]-basis or global-basis 
type. Of those of global*basis type we mention the various versions of the 
spectral method and the p version of the finite*element method. Finite* 
element methods, considered in more detail below, and coHocatipn methods 
are typically of the local^basis tjiJe. We also note that variatiopal methods 
can be adaptive. 




4.1*3 Finite-Element Methods 

Finit^e-element methods arise if for S and V we choose spaces Sh &nd Vh 
of piecewise polynomials of fixed degree Uhat is to say, the underlying 
spatial domain for the problem is broken up into small geometric pieces, 
called "elements,"^ whose size is measured by a parameter The functions 
ui^ both Sh sJ^d Vh are then restricted to be poljmomials in x and y on 
each piece but allowed to be different polynomials on the di^erent pieces. 
In designing finite-element methods, i.e., in selecting one 
attempts to achieve approximability, stability, and ^sterns of equations 
that can be solved effectively. 

Approximability here refers to the ability of the ^ce Sh ^ ^p- 
proximate the unknown solution u. The solution u is unknown a priorij 
and often only the information u E H is available. In such situations 
^ Sh has to be selected so that every function in H can be approximated 
sufficiently w^ll-by one in Sh- However, an approximation based on a few 
large elements can provide additional toformation on % which can be used 
in turn to refine those elements. This type of gdaptive element selection 
is especially important for problems with sharj^y varying solutions (see 
Section 4.2). ■ ^ 

In the choice of the bilinear form B one has, in effect, the freedom to 
choose a variety of variational principles^ many of which have a natural 
p^iysical connection with the origmal problem. For instance^ the stan* 
dard Ritz method is based on the'principie of minimum potential energy. 
An alternate variational principle is the complementary energy principle. 
Apprqpdmatipn methods ba^d on this principle are referred to as com* 
plementaiy energy or eqijilibribum methods. These methqdijrtVolve a 
constraint that is usually difficult to satisfy. One way to dr^pnvent. this 
difficulty is to treat the constraint by means of a^|range multiplier* This 



leads to the so-called mixed methods-^ They appear to be promising for 
many important problems and have recently received a large' amount of 
attention* 



4.L4'IVan9formation Methods 

In the transformation or pseudospectral method, l;lie discrete approxima- 
tion u is fint mapped by a transfomialioh of the form 

■ j=i 

into the -j^-dimensional function space of the coefficients, ay. The basis 
functions and the transformation are chosen so that T and T'^ are 
fast (order of m log m operations) and so that differentiation D is simple 
in transform space. The derivative a^roximation can then be written as 

Some common transforms are baaed on the fast Fourier transform, >vhere 
the 0j are trigonometric functions or Chebyshev or Legendre polynomials^ 
Selecting the 0j as piec§3vise poljmomiats vdth compact support> such as 
the B-9^1ines> is another good choice. By choosing the transfonnation 
to incorporate some crucial property such as the periodicity or symmetry 
of u one can improve the accuracy^ of the method for a fixed number of 
basis fuflctions. This can sometinres besi be done by choosing the basis 
functions close to the eigenfunctions of the differential equation. 

\ ' ^"^ 
4<1.5 Semidiscrete Methods \ 



When solving a partial differential equ^tioa<^9 sometimes discretis^es ^th 
respect to some but not all of the variabl^. For 6xample» for the diffusion 
equation governing the cooling of a hot rodume may discretize vdth respect 
to, the space but not the time variable^ i^ereby replacing the original 
partiad differential equation by a ^stem of ordinary diiferential ^quation^ 
Such an approximation method is r^erred to as^ semidiscrete method or 
as the method of lines. Semidiscrete methods maj; b^^used for.hyperbplic 
as lyeli as pai^abolic equations. 



One has, of course, eventually to discrettze with respect to the time 
variable as well, Semidiscrete methods are based essentially on the view 
that very effective time discretization methods are available (in the form of 
preprogrammed software packages for the solution of ordinary differential 
equations) and that the spatial discretization is the main concern. An 
alternative point of view is to consider both discretizations simultaneously. 
Such fully discrete methods have been analyzed and tested extensively. 

Among important research problems for semidiscrete methods we 
mention the problem of adaptive mesh selection for the spatial discretiza<- 
tion. 



4.L6 Method of Characteristics^ 

This is a method for hyper^oliq> equations, particularly for those involving 
only one space dimension. In these equations the solution at some point 
in space-time depends primarily on its values along a fixed, finite number 
of curves (characteristics) going back in time from tbe given point. The 
approximations to t^ese values are determined from difference equations 
that are closely related to the characteristics of the differential equations. 
This method has a' natural generalization to quasi^Hnear second^order 
equations in two independent variables. It is especially important for 
problems whose solutions have shocks. 



4,2 ADAPTIVE GFIID METHODS 

For realistic problems it is rarely feasible to design numerical processes that^ 
are reliable and accurate at a reasonable cost, and yet that do not utilize 
some form of adaptivity. Put simply, most two-dimension^ and virtually 
all three-dimensional problems are undercomputed without this technique. 
This statement vill abnost certainly remain true after the next generation 
of computen is availably The adaptive approach is to distribute the com- 
putational effort nonuniformly, so as to concentrate on the most singular 
parts of the solution or the most important parts of the solption^rom the ^ 
point of view of critical design parameters, Correspondiiii^y, ^ne must 
give less^computational effort to the regular or less critical parts of the 
solution (also see Section 4,9). At the same time, adaptive approaches^ 
may also lead to a simplification of the input data needed for the program 
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and henc€ free the user of part of the drudgery typical in preparing such 
input aiid in having to make the many a priori decisions required by most 
of today's, programs. 

The goal of adaptive grid methods is to choose .a grid that is par- 
ticularly refined, or that is aligned or oriented optimally^ with respect to 
the solution in regions of space and time that^^ critical for solution ac- 
curacy. Thus adaptive grids utilize local mesh refitment or optimal local 
mesh orientation. 

The simplest adaptive grid is a preliminary, analytically determined 
coordinate transformation. For example, in the transformed coordinates 
the problem may be approsdmately independent of one variable or other- 
WUe> simplified. The next strategy is to choose a numerically deter- 
mined coordinate transfprmation. Topically, in two dimensions the coor^ 
dinate transformation is obtained by the solution of a subsidiaiy partial ' 
difTerential el^uation. The resulting grids may be expected to give both 
improved mesh refinement and mesh orientation. The method ia^some- 
what problem dependent and can give rise to discretization errors in the 
mapping of solution variables between the various grids. 

Lagrangian grids for time-dependent problems are adaptive for 
material interface singularitieSp because the grid points move with the 
material particles. Since these well-established methods also develop 
^ rezoning and mesh entanglement problems, they provide a reservoir of ex- 
perience concerning the difficulties assodated with other evolvi ng, adaptive 
algorithms. . 

A refined grid can be constructed without recourse to a coordinate 
transfonnation. In response to a critical solution parSbieter or solution 
error criteria, selected regions of space can be flagged^ preferred orienta- 
tions selected, and refined subgrids introduced locally. Then^ small time 
steps are chosen on the finer subgrid, and an interpolation problem must 
be solved to blend fine and coarse grid solution values. Finally^ the con- 
struction is recursive, so that three^ four^ etc. levels of refinement can be 
introduced automatically^ in response to som& a posteriori error estimate 
'on the^nekt coarser grid. Predsely defined, reliable error estimators appear 
to be central to the design of effective adaptive procedures. 

In the context of finite elements for elliptic problems there is a devel-^ 
oped theory. for a posteriori mot estimates (axrdj hence, for^aptivity), 
which is based on local analysis. For a given mesh made up of elements 

...f An and with the aid of the corresponding approximate solution^ an 
error indicator rj^ is computed for each element A/. For certain classes of 
problems it has been proven that asymptotical^ the errors be$^me optimal 
when alhthe i^j become equal. This equilibration principle providers the 
basis for^he control of the adaptive process. 'Iif essence only those elements 
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Aj are refined for which the error indicators are too large in comparison 
with those of the other elements. And elements thai are unnecessarily 
small are combined with their neighbors. The study of suitable algorithms 
for this has only recently begun. For example, some results have been ob- 
tained for algorithms based on the assumption that during the refinement 
process none of the elements could be combined without increasing the 
maximum value of the error indicators, ' ^ > 

In time^dependent problems, differential equations can be ihlroduced 
for the ever-changing optimal location of the grid nodes. These equations 
are then solved simultaneously with the original parti^ differential equa- 
tion, leading to what are known as moving mesh methods. 

Alignment of the grid without refinement is also possible, if fixed 
or moving boundaries or interior interfaces are specified as part of the 
problem. By alignment (without refinement),' a regular grid index structure 
c^n be^preserved. Maintaining the grid structure has the advantage that it 
potentially allows fast iterative methods, such as fast Fourier transforms, 
to be used as^part of the solution algorithm. 

Adaptive grid techniques have been applied successfully to a range of 
time*independent problems, .The newer time*dependent methods need to 
be developed to the point where they j^an be applied, to meaningful two- 
and especially three-dimensional problems*^.d compared with alternative 
methods. An important question requiring further attention is the con* 
structjon of efficient, computable a posteriori error estimates for realistic 
classes of problems. Specifically, even for steady-state problems, there is 
the question of the validity of the equilibration principle mentioned above 
as well as the design and aiialysis of suitable adaptive control laws to 
implement this principle,^ This latter problem may require examination 
and incorporation of results in such fidds as automatic conUt»l theory, 
artificial intelligence^ and learning pn>cesses. It also raises the problem 
of the choice of appropriate refinement techniques that' produce medhes 
with desirable properties. Last but not least there is the question of data 
management, yfhich must be handled successfully to control the vastly in- , 
creased algorithmic complexity associated with adaptivity aaid to achieve 
computational efficiency. 
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4.3 METHODS FOR SOLVING DISCRETI^ DIFFERENTIAL 
EQUATIONS 

Pinite^diGTerence atid Amte^element discretizations of partial differential 
equations usually .gjve rise to laige systems of equations in which each 

. unknown is coupled to only a few of the other unknowns. Systems with 
tens or hundreds of thousands of equations are relatively comi^on. Fbr 
time^dependent problems these ^stems arise, ftom the use of implicit 
time discretizations. Fbr sufficiently fine grids, the numerical solution of 
these systems consumes a major part of the computer time for an entire 

, siipulation. ^ A 

Most nonlinear ^stems are solved hy some form of iterate method 
based on linearization^ such as Newton^s method. At each step^ these 
methods result in large sparse linear systems. In man^ cases, iterative 
methods converge only if the initial guess is sufficient cl6se to the solu- 
tiorii . 

^ There are ivfo basic approaches to solving taige sparse linearsy^tems 
of equations: direc^t/^parse matrix methods (i.e., some form of Gaussian 
elimination that vtakes advantage of the location of most of the zeros) 
and iterative methods, where the sparsily makes each iteration relatively 
inexpensive. No method is best for all problems. Fbr many prol)lems 
a combination of methods is attractive. Usually this takes the foi^ of a . 
block iterative scheme using a direct method to solve the suVsystems whose 
diagonal blocks are matrices. Combination m^thod^ aVe of particular 
importance for the solution of linear systems aris^^^from coupled systems 
of partial differential equations. ' j 

Direct methods are relatively well undersj^od today, and a number 
of excellent implementations of them exist for serial 'computers. Direct 
methods vary in the extent to which th^ take zeros into account. The 
simplest, nontrivial^ approach is band elimination^ which takes account 
only of those zeros occurring outside the band determined by the extreme 
nonzeros <i the matrix. The most complex approach is general sparse 
elimination in which all the entries that remain zero during the ^Umina* 
tion'are taken into account. Fbr systems that do not require pivoting 
(etg., those havillg symmetric, pjoaitive d^ite matrices) there is a ^rm- 
boUc preprocessing phase in whieh the location of th^e^zerq^tries is 

^ calculated. There 'exist good techniques (the nested dilsection ordering 
and the nj^imum degree ordering).for atrangfug the uijcpovms and'the^ 

' equatioiis so as to minimize the zero fill-in during elifhinatfon^ Fbr model 
problems it has been shown tl^t the nested dissection ordering pi^yides 
asymptotically optimal results with respect tp work and dtorag^^ The 



r minimum degree algorithm is a valuable heuristic approach that is com- 
petitive in practice but that has defied analysis. For systems that require 
pivbting for numerical stability one cannot compute the zeu> structure a 
priori Moreover, the ordering of the unknowns and equations to minimize 
aero fill-in will usually be significantly altered, . ^ 

^Sonie of the strong points of direct methods are as follows: (1) They 
provide an exact answer (modulo round-off error) to the linear system with 
a fixed number of operations independent'^of the ^stem's condition num- 
ber. Most production structural analysis packages such as NASTRAN use 
some form of direct method^ even for three-dimensional problems. In these 
applications (many of which are based fourth-order elliptic problems) 
it is necessary to use higher-order precision because of conditioning prob- 
lems. (2) For problems with two dependent ^atial variably, their execu- 
tion time is often shorter than the execution time for iterative methods^ 
especially for problemB of ^moderate siae^ Some of their principal disad- 
vantages are the following: fl) they require 03nsiderabiy more storage than 
iterative methods^ even for two dimensions; (2) they will almost always be 
noncompetitive with iterative methods foj^ thr^e^dimensional problems in 
terms of running time;' ifiid (3) except for the simpler'methods such as 
b^d diminution (and to a lesser extent profile or envelope elimination) 
they do not vectorize well. In fact, owing to^lhe necessary indirect aid- 
dressing involved^ there are as yet no efficient implementations of genial 
sparse direct matrix methods for the CRAY-1 or CDC CYBER-205 super- 
computers* 

Except for structural problems, iterative methods are commonly used. 
Classical iterative methods^ such as successive overrelaxation^ Chebyshev 
semi^iterative methods^ and such nev^r methods as the preconditioned 
conjugate gi^adxent method, are fairly well understood for ^mmetric, posi- 
tive defirilte systems^ and they are easy to implement,^ However^ the situa- 
tion is not so bright for nonsymmetric or indefinite systems^ though in 
practice classical iterative niethods msy^ converge rapidly with a clever 
choice of the iteration parameters. Nonsymmetnc. systems with indefinite 
symmetric parts are espepially dijfiicult to solve, and no satisfactory general 
algorithms are l^nown.at this time. Such problems arise in the simulation 
of secondary and tertiary thermal recovery techniques for petroleuni reser* 
voirs {see Section 2.^). . , . ^ ' 

^ Some of the principal advantages of iterative me^ods are, the follow* 

U They tend to retjuire minimal storage (proportional to the number' 
of uiiknowns): ^ ' - - / , / 



2. They are reaaonably fast for a vnde range of problems. Moreover, 
the number of iterations requi^d is independent of the number of space 
dimensions in the underlying p&Hial differential equation (t^ut not of their 
oii^eC] or of the mesh size); 

3. They can take advantage of good initial guesses .to the solution, as 
would be available in timfrdep«ident or nonlinear problems; and 

4. Many of than vectorize reasonably well on supercomputers. 

Some of their disadvantages follow: ^ * 

1. Mathematically rigorous stopping criteria may be difficult to for- 
mulate, e.g.,.for linear, systems with matrices that are not symmetric, 
positive definite; 

2. Many of the methods require a selection of iteration parameters, 
and the performance of the methods depends critically on such a selection 
(this difficulty is b^g overcome somewhat by the rdatively new adaptive 
methods and the preconditioned coiyug^.gradient*type metbods}i 

3. The interaction between linearization and iteration is not well 
undlrstood, especially for discretizations of coupled systems of partial 
difierential equations; and 

4. Nonsymmetric or indefinite problems may cause great difficulties 
for iterative methods.. * 

The r^tivdy new multi^rid iterative method combines the well- 
understood behavior of a given iterative technique with the fact that an 
underlying differential equation is being solved approximately. It alter-\ 
nates iterating,, toward the differential .equation's solution, on fine and 
coarser sul^ds of the spatial network, with the goal of performing no 

. more computational work on the finer (hence expensLve).gnds than is ab^ 
spljptiely nece^^. In„inany esses of practical interest, such as in jieutron 

^ diffusion in complex environments, .the technique yields sufficiently ac* 
cifrate solutions to the equations, ia^a compijlitational time proportions! to ' 
the number of unktiowns.^ This has'been a long^sought^for go^sl in ihe 
approximatioaof elliptic equations in more than one space dimension. As 
an itierfktive methqdj^.also has a natural extension to nonlinear equations^ 
and its logical structure^ together, with.the already necessary calculation 
of residual errors, points toward the ihcorporatiob of adaptivity for the 
nesting spatial ^ds. ■ * 



4.4 CONTINUATION AND HOMOTOPY METHODS 



For more than a century .the so-called continuation principle has proved 
to be aD important tool in niatheinattcal theory. For example^ early uses 
date back at least 100 years to H. A, Schwai^^and his work on conronnal 
mapping?; then, in the early part of this c^&ry, J. Hadamard and M, 
Levy applied these techniques in connection with the inversion ofnoiilmear 
mappings, and it is also a basic tool in X Leray.and J, Schauder's work on 
difierential equations. But it was essential]^. only in the early 19503, with 
the advent of automatic computing^ that continuation approves began 
to be used in numeiici^ applications, Since then they have grown into^ 
extremely powerful technique for the niunerical solution of wide classes of 
nonlinear problems, ^ 

One of the problems to which continuation techniques^ applied 
concerns the solution of a nonlinear equation F{x) ^ 0 in some apace 
X. in order to compute a specific solution z*€X a possible approach is 
to imbed the equation into a famify of equations H(x,t) = 0, 0<tpl, 
for which there exists a computabfe solution path x = x{t)f 0<£< f in 
X beginning at a known point.^0) and ending at the desired solution 
x{l) ^ x^. in other words^ one considers a family of smoothly changmg 
problems, the final problem being the oi^al problem in q[Uestton and the 
initial problem being one whose solution ia easily determine_d, Use of the 
continuation, therefore, requires an ability to solv^sequence of problems 
wheirthe solutions to n^by problems are known, 

A related^ but conceptually distinct^ problem ari^ in maiiy applieft^ 
,^ions in science and engineering where the equation under consideration 
always includes a number of physically relevant^ intrinfiiC parameters, in 
other words, the equation has the gen^ form F(Xfp) = 0, where x belongs 
to some state sp^ X and p varies in a parameter space P. in this setting 
it is rarely of interest.to^ejiermine only a few apecifit^ solutions. Instead, 
the requirement is to follow paths on the solution manifold € Xx 

P; F{x, p) = 0} in the space XxP of all state and parameter variables and 
thereby to detect specific features of the manifold that signify,a change of 
behavior of the systemjinder study,. 

For the first of these problems, namely the computation of a specific 
solution oT some question, two distinct continuation techniques aie^^vail- 
able,',The first involves the case when the path of solutions to the family 
is smooth, which in turn allows its representation as a^ solution of a 
difTerenttal equation. The second approach is based on liomologieai rather 
than homotopic concepts and makes use of,a numerical form of 3 lesult in 
algebraic topology (Spemer's lemma). This approach has been reformu- 
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lated and is now considered principally in the formof alg)c>rithiiisinvej^Qg 
piecewiee^linear solution paths. MucH of the research in thl^^area 
concerns iiiese piecewiee-lmear con^muation methodland their appllcatto;i 
t^fi^ed'point problems in economics optimization and game theory. 

Continuation methods for foUowing patbson the solution- manifold of 
parametrized equations developed mainly m structural engineering under 
the name of ificremental methods. Applications to other areas, as, for ^ 
exampfe, to computational fluid, dynamics and to phase transitions in 
statistical, phjt^sics have only begun to appear reiativjely recently. For a 
numerical analysjs of ^given solution, manifold continuation methods h^ve 
to be considered^in a broader sense as a collection of numerical procedures 
for completing a variety of tasks, including the following: 



^ _ fbJlcFsy num^<*ally any^urre on the manifold specified by an,o r f"^ 
^prwri j^ven combination of parameter values with one degi:!eeof freed9m.l^,^_: 
\ , ^. 0^ &Ay such (turw debernjine the 6xact location of target points" " 
.'^^berc'^J^vei-stete- variable valiie* / .r" , 

J ^^^TiieS&^jatOTKl^ 
■ sdliitloilt^ja^pji - ^'}y^\[:Z' ■ V ■ ' ' ' 

4." Oir^^y curve' idenl^y and eompute\tbe initi^ points where 
■_ Ability may'be-Jtfst;^ , - - . , " / - - ^ > 

.5.. '^^jn^aiiy pne of, the critical points foflow a patli in the.critical . ^, 
boundaiy.' jln."con&ast to the cas^. under % such j^ths .are t55)ic^ly \ 
sp^cdfiedtjbjr Cdi|^^bi^^ pffl;amejie^^th:t^ degrei&s of Xreedgm 

' tc^ti^ mihj points of.th^ pathare,in 

£h^.«ntii^^ ! * ^ 

; &. On aiqr>plution path ^e^erminethe ItySion of biTurcatignTppints . ^ . i 
' ,and^"tbtp^ths.intera,^tirig-:^i^iIJ^^ /--^ J. - / . ^ -^^C^ i"- ^ . ^ ' -^^.j^^' 

I \" * -'"V ""^i/i^ \ ' V V \\ . - * * 

F6t Specific apjpiicaiipns adctftionjal taslca may aria^., i^stan.ce, iO, v^^:^^;^^ 
the p^r^etrizeci equation repi^^sQits th^eqp^^ , , . ^ ^^^--^ 

ofdifTeitittiEdf^ationsthenW^^ay^^^ . .-^^ 

- on a particular. solution j^ath |>e]^5^ . , ''-^ - 

;\^y^ml)ranch off"^ stj^tV^f^^ i.^' ^ ■ :'" 

In recent years much res^^<^ has bn^n <t^^d to ihese v^us*^^ . " 
bt& ther^^r^ain many open jqpestion? eaft^siaJfe: in^jyJimectiop. vri^tl^^^^^ , _ 
mor« spcclaiiied Ja^^ CTLcK?^ JheJ<%t!^j^^pf ^ ^ 
^ Mor^ver, Ubraiy-^uali^.p^^ "!i 
Sevelopmeijtp ahjlthc w 

pf ^qRv^re fer tjig solution pf or . ,,V- ^ ^T, 
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"apptications, as, for examplj&j^ to fluid-dynamic problems or to the case 
of-Structutal probl^m^jiiyolvmg plasticityt creep^ and yiscoelastic effects^ 
the situatioii is stilL^id^'opea^ and only relatively ad hoc techniques are 
.■4>raaJabIe^".^_ , * .^i. 



, , - / 4.5 OPTIMIZATION METHODS - o 

'J'iV ■■ ' / * ' 

, r / Opftinizfiltion problems occur in all areas of. fffiieiicer engineering, 
IScs^^} ■ ^economics, and applied statistics^ They^may mvolve some .least-squares 
.appi^)xijn.aMQQ <3{jDbse occurring m a math- 

^^V .^,. - esoa^icd model^oh the basis of 3e%peiim 

1, / - the desigiT^^a^engineering^structure^ optimal control of an engine^ng^ 

... ./ , ^. .systejp^^ori^^^eling ,of economic or busm^^ystems.. These are only a 
___^.^-j,. -/ew exam 

^ Bniadljr;§pcakii^ aU^^pf these pipblems a reaJ function, ^^pally _ 
. . _ , . ^c^iiedjiije objective I^Ettietipn, to be miiumized or maximized over some 

i . ; ^^/^jpqgst^^ntset'in a The problems 

on the mathematical charactcris^<^s ^^tfefe.^ 
>-.i>C7^;> v^^,^?|^fe constraints^ the dimension of the ^pac^*TEe ^. 

'Z^ . ^Jj^l^" "a^tint^irf^cojmp^i^ that is available, ond the requirements 

^ " A:"^oa:0^^mew, of the research problems jt^ncemmg computa: „ - 

.t^nal met hods for such optim^aUonprobJema in fiiute-Sm^ 

i^ien^^n ^ r^he rfpqrt Program Directhnh for ^Computational ^^-'V^ 

' * ~"Yi^^T^ia^^^ Department, of Enei^ (edited ^by fi. 




T^i^. — £_rL-:;_..r_--t.^ (JeaJing with nonhn^ 



.- T^.r ^.r- Production oChig!l:;9^ity software and-other software-directei 




5. Techmfliififa'tol«s^^^^ /* 

7a 'Metfeods"I6r:non3iffei:e^ 
separable or con§;ttftili^ least ]3<Bia^,^^^ 

altematiVniethods"for,fi$4Wng.nQQ^^ / 

* 10. Further study , of ,3^e{lgd |Kp©£t& ^IK^^l^^ 
"such as special *yp^,ttf^ieonsW^^^ 
techniques, the merits of vmj>y^ :ft^iye . set, \- 

or scaling. " V:A 

■ . \^,> : .V. ,^ 

T^ie report does not address lnfini^:4itn^<e^ 

mathe^madcAl ptpblemj o*^. l>€ca]^|pria^ ^^^^^L teymiques ^ - \. 

areln^edrDmostappiiQatiQji^/F^^^ . \ . 

6ipal approach ^.^^-l., 
cretized form .of Uij^ jfir^^^ 




Many phof^cai ^^^^ . 
^ 5^ diWr^te_d|^^^ 

,^ coi^ipdSt^e5.,aiid.iifl^i|n^^^ . _ 

^ t^W fijii 

The m^dci^i^Dlf B^^^ " 
, JsjiaUier natu^Sl/fpr^^^pT^j.lft^^ 

" Inie rb^t^jll^^^^^^^ - 




Corrcsjpoxitiiiig to the wide range of applications, the fornis of graph 
problems differ widely. For certain network models interest may center 
on the connectivity pro))erties of the grapiii that ia, on the.deWmmat^^^ . 
of whether a particular com.ihp*dity, can tiijEtiispofted ^i^tween specific 
vertii^ Then questions arU^ a^ou^ ih^ maxj^^ possible flow.ihat can 
be accommodated under partjcubr constraints. ThesteLprobtemsin tum^are 
related to theso-calledViiiyiera^ity gnd reliabilfty. problems for networks. 
On the other.han^ if tim^_€nt^rV into consideration then questions of 
waiting time and of b^st rQUting].may,Mi3C. XheJse a^e onfy! ^ibw of many 

types of such prob!eriS&.77^^'\^>'\v^^^ ^ ' , 

FYom a computational viewpoint theee v^(»is problems, have stirnu-^ ^ 
^ lated the development of numerous combinaWrial and grapl^^Keoretka!^ 
goritWs. But there Femain.oiahy.opeti resc^n^Lquestbns, eSpeciallVwh^ 
it com,e^ to the production brg^enera[,^tpbust software and tlie availabi^lity^ 
of aigorithqi^^Jor pK>blemsJnvblving'l*tfgp ^ 
Graph^heD^.^tica] formulations ^ a^so finding incre^ing appiicatiqn 
in connection with, the niimerical solutloij of problems that do not have aff 
inherent graph-structtire. JCTne^such class ^of, problems concerns the com- 
putational tian(mng of large sparse -matrices^ ^any of the algorithms used 
here perform a sequence of steps tW transfonp.t^einatrix into eucces- / 
sively simpler matrices of^say, a more nearly^ up}>€r^i^ngular or diagonal 
type. These transformations achieve their aiftx by. introctucing aeros in 
place of originally nonzero matrix eJements. But, as an unavoidable. by- 
product tbey ^30 introduce nonzero ele|pents3n^p1aip;es where the origmai " 
matrix ^elements were zero. Thus^ J|)kO!t[er 'to exd^ the sparsity siruc: 
ture of .the matrix one must provide a^Sj^^^UTicluK th^ either. all&Cates^^^ 
from the o^ife^t space for all 4r^e ^J:in "dt^^ or that 

.allows for a d):&aitiical allocatioKjf sps^^e^for^ itjocci^Jn.^ 
both cases^ graj^b^tliegrettcal ap|rrqacl][^form ji^^^f^lor the design of- 
desirable alg^%h}^. *' ''-^^ , ■ V"^ ' ' "-V^^ 

An example jof ttjg^ use ©f^a static data structure is'thg.c^S^^^Tti^Ka jo)y 



and^timn peprese|iiaU(^n .can1>nngthe matrik into.}^^ 
small bandwiilth. TMsr^sTrequently the casi^T)( tjie^ 
in the ffnite-element disci^Ucatipn of elliptic boundary- vaiu^^ problen^^^ 
and the corresponding bandwidth optimlzatjoii routines are wid^ used. 
in practireyOa thepther band^if a dynamic^^torage strycturelsused^ th^ir. 
special ca^e had to be taken in the d^esigti of pivoting strategies fqrf^ucjqg 
fill-in while at the same time maintaining numerical stability. , Various 
sparse matrix paclcages hiye been developed for this purpose Th^ are 
typfially rather la?^ and cpmpl^-piec^e^ of Software. But the field is stUl.J, 
tMber active djevelopmei^^t, anc|.t^e5e.re5iaiiifloany open research que8tipjis,.J^^ 






In particular, n^cU still jieede to be done on the interrelation between 
sparse matrix algorithms and special-purpose computer architectures. 

In recent years graph'^theoretical approaches have also found increas^ 
ing application in the design of algorithms for the numerical solution of 
problems of mathematical physics, in particular, of fluid dynamics. For ex- 
ample, it has been observed that certain natural finite-difi'erence discretizap 
tjons of th€ equations of viscous^ incompressible flow admit interpretations 
as systems defining flows on certain associated networks. Typically, the 
network nodes correspond to the idealized control volumes represented by 
the mesh points of the finite-difTerence equations, whereas the network 
arcs correspond to the paths on which one may identify the discrete finite- 
difference mass fluxes. Such observations can lead to remarkable savings 
in computing costs because they open the way to a priori transformations 
.of the originsJ (di^retized) system of equations to completely equivalent 
systems in substantially fewer variables. For flow problems in two and 
, three space dimensions the^reduction factors are nominally 1/3 and 1/2, 
respepilvely! In graph-theoretical terms this approach corresponds to the 
.'^transformation Qf the Maxwell-node equations to the Maxw€ll-mesh equa- 
, t^ons, fong used in electrical-circuit problems. . . 
^ ^ The jkey rejM^ii^e^ 

^ flow problems is that the discretized momentum and continuity equations 

interpr^table as "Ohm^s laws" and ''Kirchhoff-node laws" on an as^ , 
;^ ,s^wajt^ This requirement p^mtts great generality in the actual 

^^ .y^ ^ fono law3 and easily acf:bmmodates both implicit time- dependent 

/ As w^lVad nonlinear steady-state difference forms of the Navier-Stokes • 
;:\ \€0Uiattonsy^ > , , * * ' " 

] V ^|*3^y^sPp?sib\^.*^*^ 
\^ i^^.Py?^^%^?;^ H^J^.^jh^new idea aH)ear3 to be the introduction , of pseudo- 
^ !,!Jo^9 iptp ^^Jje tio^e {itw.and the identification of corresponding pseudtf* 
V^^rAel^^^ the Ohm *s lawa obtained ffom, the momentum 

-'^^1 felpation^,^^ ca^sej^^this prpdu<jes numerical schemes \ 

irtde^iSentof Ipem^cQustic velocities. 
^^^'^,^^^^T^e jep^ development. For ex- 

^l|,^^?i^e,^th^^^ prob^ns und^^the 

..^tV^^Wryiej^jS^^ Important special case, 

^^^^^S^^)^' P^Spl^^^^ Jthe design of reliable 

^^ '^^^li^^^^^^yig^sjwte^s tJiersaJ^y aniiysi^ of nuclear power plants. 
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4J SPLITTING METHODS AND DEFECT CORRECTIONS 



SpUtting is a means to separate a large intractable problem into a sequence 
of smaller, or at least more easily solved, problems. These methods are 
invoked to reduce significaDtly the computational efTort (time and memory) 
needed to solve a problem compared with solvLiig it directly. Often a 
form oT splitting IS required to make solution practical, and some cases 
even tractable. Various forms of splitting are common in engineering and 
scientific applications, even if they are not always recognizeJas such. 

Defect correction is a widely used, if often unlabeled technique of 
splitting- It presumes that one wants to solve a given hard problem, that 
one has in hand a guess at its solution, and that one also has a nearby 
problem that can be solved easily. It corrects the guess by solving the 
easier problem with special data computed from the guess'. 

To illustrate the method, suppose that after discretization it is nece^^ 
sary to solve a finite-dimensional system of equations of the form 

' • . • A{v)-b^(i ^ , . 

where >1 is a nofiHnear operator, is a known vector, and v is the solution. 
Often V is difficult to obtain directly, but the fesidual error 

for an approximate solution w is easy to evaluate. If there is a related 
system^ 

that approximates the original system and that is easier to solve, the defect 
correction algorithm may be appropriate. Given an estimate near a 
root Vn^i of the original system, we can expand the original equation in 
a Taylor series to get " ^ 

O^A{Vr,)'b+P(Vn^i]■-P(Vn)-^^^-J^k){Vn^^ 

where t - Vn-^\ - The defect correction itemtion is any approximation 
to the. above equation. The simplest such iteration is 

This iteration will converge* it -t^^ and (the Jacobian of ^) are^near 
enough to Vn-n and:^;^, r^'peSjtively. 



One of the most common splitting algorithms in engineering applica- 
tions is the alternating direction implicitj^^I) method- Here, the model 
problem is a partial diETerential equation system having two or more in- 
dependent spatial coordinates. Its direct solution requires too much com- 
puter storage and time to be tractable. Using an ADI splitting, the prob- 
lem is reduced to a manageable sequence of one-dimensional problems. 
In terms of the defect corrections, the system P is the sequence of one- 
dimensional problems that is much 'feasier to solve, yet that closely ap- 
proximates the original system- 

Other forms of splitting also ^se- For examf^ie, in certain systems of 
partial differential equations some of the equations are weakly coupled to 
the others. In these cases/ solving the equations sequentially (rather than 
coupled) can result in significant savings. Similarly, in some combustion 
problems, considerable savings are realized through operator splitting al^ 
gorithms in v/^ch the chemical rate terms are handled separately from the 
sp^ies transport terms. These methods are equivalent to matrix splittings 
of the Jacoblan of the systwn. 

Analysis splitting, methods important because th^ methods often 
do not converge/We must be concerned*about accuracy and convergence 
of the factored system.^ Even though each iteration may be fast, we 
must have some idea about the overall convergence rates and about any 
degradation in accuracy resulting from the split. Analysitwill also likely 
lead to^ the identification of matrix properties that suggest, a beneficial 
splitting that would not be apparent fr om pbysjcal r easoning. Splitting 
can be considered as either a splitting of the original equations or as an 
approximate factorization ofthe iteration matrix. The former is the most 
intuiUve, and the one in which physical insight is^valuable. However, the 
latter is probably the one mdih amenable to analysis. 

4.8 MONTE CARXX) TECHNIQUES 

Mathematical solution methods can he broken into two broad classes, 
Monte Carlo methods and deterministic methods, depending on whether 
chance events are included in the method. All scientists are familiar 
with deterministic methods hut most have little familianty with Monte 
jQario. methods. Deterministic .methods do not depend on chance, siid 
Calculations 'Performed using the s^e input data ^^hd the same computer 
. code will prj^de exactly tiie same results. Monte Carlo calculation^, using 
the same ^nput data ahd the same computer code will provide dlETerent 
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results, depending on what chance events occur. In this section we discuss 
what Monte Cario methods are and how they might be improved, 

. The Monte Carlo method estimates averages on a probability model. 
Any quantity that can be expressed as an average value can be calculated 
by Monte Carlo techniques. Sometimes a probability model is obvious 
from the problem'itself. For example, the probability that tossing three 
fair coins and obtaining a desired outcome (for example two heads (H) and 
one tail (T)] is easily estimated by Monte Carlo methods. The probability 
model consists of assigning a probability of ^ to an H and a probability 
of ^ to a T on each toss and assigning a score of I to a desired outcome 
(PIHT, HTH, TPIH) and a score of 0 to any other outcome (HHH, TTT, 
THT, HTT, TTH). The computer has a random number generator that 
generates random numbers uniformly on the unit interval (0^1). A uniform 
distribution means that any random number ^ is equally likely to have any 
value between 0 and L Thus a coin toss can be simulated on a computer 
by. 

H occurs if ^ < 
T occurs if ^ > 

To toss three coins the computer uses three random numbers ^i, ^2t 
and A typicW set of tosses .might be =^ OJ (T), ^2 = 0,1 (H), and 
^3 = 0,4 (H), scoring 1, The probability of ob taining two heads and one 
tail is (approximately) the average scor^l[|) after many^ets oTcomputer 
tosses* , ' * , ' 

Sometimes a probability, model is not immediately apparent, but after 
a little thought the desired calculation can often be expressed as the 
estimation of an average value. For example, the integral 

D — aJa 

can be thou^t of as the average value of f{x) on the intervaL{a,6), To 
see this, note that by definition^ the average value of /(x) with respect to 
a probability density p{z) is 

Thus / is the average value of f{x) with respect to the probability density 
p{x) ~ 1/(6 -a). Here p(x) is a uniform probability density on (a, 6), and 




X values can be samples uniformly on (a,&) by setting 



x, = a + (&-o)Ci 

whereupon 



where / is the estimate of /. 

There are many applications for the Monte Carlo technique, but the 
application to neutron and photon transport has probably consumed more 
computer time than all other applications combined. Because of this, and 
because neutron and photon transport problems have a natural probability 
model, the remaining discussion will pertain to neutron or photon trans- 
'port problems. PUrthennore> because neutron and photon transport are 
^ similar, the remaining discussion will refer to particle transports 

Particle flow jn nature is a random process. Qf^ cannot se^ exactly 
what will happen to an individual particle. One can only say what the 
probabilities are that various events occur. For example^ a particle trav^- 
ing through matter has a certain probability to collide [)er unit distance; 
the actual distance between collisions is unknown^ but the probability of 
traveling a distance I without collisions is known. Similarly, the nuclide 
a pa rticle will c ollide yrith is not known>^t the probability of colJSdbg 
"with tlie nuclideirknown. ' - 

The simplest Monte Carlo model for particle trans|]k)rt problems uses 
the natural probabilities that various events occur for the probability 
model in essentially the same way as tl^e coin toss example. That is, 
paj'ticles are followed fh)m event .toeventin acomputer> and the next event 
is always selected .(using the random nmnber generator) tvom a number of 
possible next events According to the natural event probabilities. This 
type of model is called the analog Monte Cjarlo model because it i& directly 
analogous to the naturally occurring transporjP 

The analog Monte Carlo model works well when a significant frac^ 
tion of the particles contributes to the estimate of the average. This is 
analogous to having a significant fraction of the particles in the physical 
situation entering a detector. Tljere'are niany problems for which the 
fraction of particles contribution (scoring) is very amall> l^ss than 10^^. 
For these problems, analog Kfonte Carlo fails because fqy> if any, of the 
particles contribute> and the statistical unce^tmt^ in theans}^ is imac*- 
ceptable. " j - \ f n 

^ Although the analog Monte Carlo model is probably the simplest 
conceptual probability model, there aire an infi^jjte number of probability 
■/ " . .' \ ■ . ^ 



models for particle transport .that will estimate the same average v&lue 
as the analog Monte, Carlo model These' other techniques are known 
as nonati&log Monte Carlo models, and they are useful because although 
th^ average yalue remains unchanged, the variance (uifcertaiaty) of the 
estimate can often be made much smaller than the variance for the anal^ 
estimate. Practically this means that problems that would be imfH)ssible_ 
to solve in days of computer, time can be solved in minutes of computer 
time. ' ' 

A nonanalog Mont^ Carlo model attempts to follow "interesting^' par- 
ticles more oft^n than noninteresting pftrticles. An "interesting" particle, 
by definition, is a particle that contributes a la^ amount to the quantity 
(or quantities) that needs to be estimated. There are many nonanalog 
techniques, and th^, all are meant to increase th4 odds that a particle 
scores (contributes). Tonsure that the average score is the same in the 
nonanalog game as in the anaiog game, the score is modified to remove 
the effect of biasing (changing) the natural odds- Thus, if a particle is 
artificially made [q times as likely to execute a given random walk, then 
the particle's score is weighted. by (that is, multiplied by) IJq^ The average 
score is thus preserved because the average score is the sum^ gyar all r^- 

, dqm walks, of the probability of a random walk and the score due to that ' 
random vralk- The trick in obtaining low- variance soluttons is to choose q\ 

\ such th^t all random walks contribute the same score, in fact the average, 
score- This then would be a zero- variance' solution- - ^ 

It is always possible for any (linear) j)arti61e transport problem to 
select q^s for each ^dom walk such that every partide^s score .is ^the 
average score; that is, a zero-variance solution- The hooker js, of course, 
that the perfect q'^ are ijot known, thus a zero-variance solution is not^ 
practical However, people have often been succes^l enough in guessing 
good q% th^t IS biasing the odds, to improve the calculationai efficiency 

Jay several orders of magnitude- This is-lypicaliy done iteratively vrtth a 
person making several short ^onte Carlo computer runs (calculations) and 
using information gained on each run to better guess.how to bias the next 
run- When the person^s guesses no Wger improve the calculation, a long 
run is done with the optirnal biasinP^amed in the short runs. 

Can the computer lea^ to optimize the biasing? - The ccn^^ter is, 
after ail, supplying the informatioa that the person uses to learn. la 
the-past decade a number of computer learning ^hnjques ha^ been 

. tested. ' Thus far it ^has proven impossible for a computer to- take an 
arbitrary transport problem and, without human aid, optimize the biasing- 

,However,.two thinjgs^are worthy of no6e. First, computer learning aided 
by human judgment appears to be sul>stantiaily better in many cases than 
human learning aione. This typically works by having tUe computer inform*^ 



-the person how the computer would bias the subsequent run and having 
the person selectively alter th^ computer's suggestions according to the 
person's intuition. Second^ the amount of ^human judgment required is 
decreasing. The day may come when the computer needs tio human aid. 

pnce .human aid H no longer needed, the compute^ ^n learn to 
adjust its own biasing particle, as the calculation proceeds. An interesting 
implication of an' adaptive Monte Carlo Schmque is that the common 
central limit theorem of statistics that would constrain the accuracy of 
the calculation to decrease as the square root jit^e number of paJrticles 
followed no longer appli^. The <;ommon centre? limit theorem applies 
only when each particle's random walk is independent of all others and 
the sampling process is identical for all particles. Consequently, the rate 
of convergence may be more mpid. Although some limited statistical 
theories exist for dependent random variables, it has not been investigated 
whether these theories profitably be applied to Monte Carlo particle 
transport problems. Tnu^^ is uncertain how fast an adaptive Monte 
'Carlo calculation is converging or indeed what the ma^dmum convergence 
rate for a good learning proems might be. Empirical evidence shows that 
the convergence can be substantially faster than t^e square i;oot of the 
number of particles. In light of this, Monte Carlo methods can be expedited 
to become more competitive with deterministic calculati^ms. ' 

4.9 PROBLEM-DEPENDENT METHODS . ' » 

A variety of adaptive methods have a c<>mmon goal: to overcome the 
computer size and grid rjesolution limitatipas, which are especially severe 
in singular problems/ by use of computational elements that fit or model 
the singularity more directly. In this way, it may be possible to incorporate 
into the solution algorithm considerable analj^tic information beyond that 
provided by the equations themselves. This theme occurs in maiy aspects 
of nun^ical analysis! In the finite^element method, on^ may choose 
elements that include any singularities in the solution being computed (see 
Section 3.3). Grid adaptivity {see Section ^.2) iS^glso problem dependent. 
,Bere we discuss examples of somewhat more special methods. Of necessity 
they^apply only to classes of problems that poss^ related singularities. . 
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4.9.1 The iriemann Problem and Nonlinear WVe Modes 



The Riemann problem for a nonbnear h> perbohc system of conservation 
lhv,s IS the Cauch) problem in one space dimension for data that are 
ever>'%here constant except for a single jiimp disconttnuit) at the onginr 
The solution of the Riemann problem provides a resolution of this discoh- 
tinuil> mto the nonlinear .modes, or leaves* which propagate coherently 
in tmie. Thjs idealized pn>btem can be thought of as 'an approximate 
description of a higher-dimensional fiov, field esp^ciallv in the neighbor* 
hood of a djscontinuit) hypersurface. This point of vie;^ has led to a num- ^ 
ber of improved finiie-difference^bemesi i^hich attempt to determine the 
various nonlinear i^ave modes at each point in space and time and to ad- 
just the djfTerehcing of the dififerential equations to take advantag? of this 
knowledge. This adaptivity is to the structure of state space* tn contrast 
to the coordinate space adaptivity^iscussed in Section 4.2. 



4.9,2 Front Tracking ■ ^ * 

• 

FYont tracking is a combination ^ adaptive grid methods tth the u$e.of 
Riemann problems. The method is adapted U> probflrhs that hav^e inipor* 
tant singularity h>i>ersurfaces (lines In two space dimensions^ surfaces in 
three space dimensions), such as shock waves and contact discontiniuties. 
In one version of this method, there is an overall time-dependent coordinate 
transformation to map the discontinuity into a fixed grid hne in sonte set 
of computational coordinates. In another version, the discontinuity is rep-« 
resented b> a moving lo;^^r>dimensionaI gnd^tbat follows ('*tracks") the 
^dynamical motion of the discontiiiuit)/ The motiocTof the discontimiity 
and of the movjng grid points that track it are governed bv solutions o£ 
Riemann problems, or equivalent!) b> the method of char&ctenstics. 



4.9.3 Vortex Metliod , ^ ' 

m 

This method introduces small lifieS^r surfaces of vTorticit) into a fluid flo;*. 
The method is adapted to the study of shear^laycr discontinuities^ bound- 
ary layers, boundaiy-layer separation, and turbulence. The equations <)f 
motion of an iifeal fluid yield simple equations for the motbn of a collec- 
tion t>T vortices imbedded in the fluid- fact, the, vortices mo\e pas5jvel> 
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with the fluid, and their mutual interaction is described by a Haniiltonian 
system of ordinary difTerential^equationSj with Coulomb typp interaction 
energy. In the case of the Navier-Stokes equations^ the vortex motion also 
contains a difiusion term. 

These methods have been applied successfully to the problem of tur- 
bulent flame propagation (see Section 2.2).. In this problem, the turbulent 
mixing is a priniEry factor in determining the flame speed. The turbulence 
comes from the boundary layer and in the boundary Isyex is calcuiated b} 
a \ortex method. Related methods have been developed under a variet} of 
names such as boundary integral methods and Green "s function methods. 



4.d.4 Scale Invariance and the Renonnaiization Group 

Scale transfonnations are the transformations x ^ sx^ acting on space or 
on space and time. A function u is homogeneous of degi^ a if 

u{3X^ sy, sz) = 3^ u{x, y, z) 

and scale invariant if a= 0. Many probtemsVave solutions that are scale 
mvariant or approximately scale invariant over some parameter range. 
Such solutions are called similarity solutions. Using the scale invariance, 
, oce of the independent variables can be eliminated^ making the solution 
mere eiementary to compute. ' " ' , ^ 

Howevert scale invariance can also indicate the occmrence of com- 
plex pheEjpmeQa, SpeciSca%» ai^ singularity that occurs in a scale^ 
iQvat;iant problem must be repeated in all length scales (for which the 
scale invariance holds). Mathematically^ C^tor sets^ snowSake curves, 
and fractals a^ examples of such pt^nomena* in statistical^hysics, criti- 
cal phenomena in the equation of ^tate is a scale-invariant problem. One 
general picture of turbulence holds that scale invariance (vorttces on a 
large rang^ of length scales) describe themertial range, oi^ energy transport 
"range, of turbulence. ' ^ , 

To implement these ideas in a f:omputational algorithm^ one integrates 
over a given set of length scales tn a standard manner. IJhe result of this 
cpmputation is taken as data for &^new calculation over a new set of length 
scales (with the original degrees of freedom eliminated from the prob- 
lem). This process is iterated and if convergen^t gives the^scale exponent^ 
a. Sample numerical calculations of this type were discussed in Section 
2.1 in connection with turbulence. The method is well established for pe^ 
turbative calculations (with a small parameter) in critical phenomena in 
statistical physics, ' * 
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FIGURE 44 Tfie evofetion ind merger of isolated vortex structures as 
predicted by contour dynamical techniques. [FVomi;* A. Overman Hand 
,N. J. 2abu3ky; Pky$- Fluids 25, 1297-1305 (1982).} 




Since the reiiormMI»3Ition group ^methcxl^ are noyel, it is worth men- 
Uoning that the mat^maticCl foundations of thi&lhethod have been estob- 
lished recentt> in several cases including examples of hierarchical models 
in btdtisticat mech^nicSf mter\al maps, and renormalized quantum fields. 



4.9.5 Contour Dynamical Methods 

ConjLour dynamical methods are being applied Jbo a variety of incorpipress- 
ibie f)o^vs in two dimensions. These generalizations of the 'Svaterbag"" 
niethod provide simplified models for following the evolution of contpurs 
separating regions of piecew^s^constant (pc) density that are the sources 
of the flow. The flow ia the result of the sejf and muttlai interaction of con- 
tours that evolve, mainly tsy area- preserving maps. These jjfiethods have 
been applied to the Euler equations^ ^vhere pc finite-ar^ea-vortex regions, 
and/or vortex sheets at/iensity interfaces are sources of the flaw; and 
the equations for a weakly ionized and strong!^ m^'etised ionospheric 
plasma ctoud, where pc ion-density regions are sources of theJow- For the 
former, a large class of 3teady*state translating and rotating solutions with 
pc vorticity (A'-states") have been founds Figure 4.1 shows the merger 
an^i-breaking for a perturbed corotating V-state (two pc finite-area vortex 
regions with the same cireulation)^ a familiar process in free-shear layer 
experiment^. Notice how the two regions merge to (orm one region (by 
snipping out the common boundary at t - 10) and then stabilize by Meet- 
ing vorticity in thin filaments. With these methods it is possible to obtain 
detailed information about the regions jDecause the climension is reduced 
%^ One. The curvature of the contour provides a predictive signature' of 
the e\ otving small-scale structure, eg., the rolf up of vorticity filaments, 
'etc. 



4.10 ( ;OMPUTER SOFTWARE A^fD HARDWARE 

As previously discussed^ large-scale scientific calculations that tax th^ 
resources of the most powej^l^ computers will continue to be essential to 
modern research and development efTorte, To obtain long-tenn reliability 
and stability of future applioations codes, implementing and testing of 
high-level numerical software should be coordinated. This will require a * 
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strong research and development effort witKoog^tion supported ^ong 
sppllcatlona progra|timers^ the theoretic^^^Mil^^cal methods researchers^ 
atid the computer-science software developer 

RepetitioiT of expensive, error-prone, and time-consuming coding of 
commonly used methods should be avoided. Much of the current scientific _ 
software nc)w being developed is redund^t. If the common elements of the 4 
existing codes were av^Iable as modules, future applications programmers 
could use these routines and eliminate -much of their efforts. New softw^e 
is most readily '^ccepted if it is compatible .with existing techiiiques and 
simple enough so that potential users can obs^e tangibl> better results 
in a <^al run than those existin|^ methods can produce. If soth^hl^h- 
level routines were available^ they could perform many of the common 
procedures found in applications codes, including ^id generation, rezon- 
ing, numerical interpolation, differentiation, a^ integration; they could 
approximate differentta^ boundary conditioiuiand solve large, sparS^ non- 
linear systems of equations, * 
^ An important goal is the machine independency of appUcations pro- 
grams, This is a difficult task because methpds tailored specifically for a^ 
partictilar machine architecture will probab|s^become more the norm than 
the exception. We can, by keeping machine-dependent codes in libraries of 
high-level software with standard user interfaces^ strive to keep the user*s 
scientific applications codes portable. The underlying^ library routines will 
be, necessarily, less portable because the architecture of the n^ machines 
will certainly be different from that of teday*s supercomputers. To uti- 
le the inherent power? of the new machines we will have to re-examine 
traditional metWds an^ identify JLhe better ones for a^di^i^ul^r machine 
architecture/ - ^ ^ ' ^ 

The continuing revolution in microelectronics is having a profound, 
impact on scientific computing. Indeed, it is likely to change our con<^pt 
dramatically of what scientific computations can and c£^not be effectively 
performed. Mo^ certainly the Impact of this revolution will be much 
greater than, say, the impact that fioating-point hardware had. Moreover, 
viMe the costs of individual tasks will be greatly reduced, the domain of 
scientific computations will be greatly exp^tstted aAd "frontier computa- 
tions" will continue to be^^xpensive. ^ 

Thege* chfinges are belng^brought about by :a number of factors^ 
Individual components are becoming increasingly faster and smallgr^ Very- 
large-scalg integration of such components is not only reducing th^ siie of 
the packaged systeins but also providing opportunities for customized in- 
formation processors. Also, the decrease in the cost and size of computer ^ 
memory implies that we can look for much larger memory aystemsr This 
wril obviate many of the existing difficulties with secondary jstorage. 
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Despite the fact that eomponents are becoming faster^, the mnits of * 
raw madiine spee^ar^ not visible^ and further gains will eventually have 
to be made f>y using clever architectures and algorithms. Some kind of 
parallelism^^ms to be unavoidable. The programming issues involved in 
using parallel m^hincs are still not satisfactorily resolved. The automatic ^ 
detect)on of parallelism and the resulting scheduling of multiple resoijrces 
are important open problem^. ■ . ^ 

Architectures, such as systolic arrays, ba^ oiii specific subtasks cafi 
increa^ie the performance of systems involving thes^ subta^ks^ by* several 
orders of magnitude and clearly have a'bright future. Algorithms contain' 
ing compute-uitensi\e subtasks that can be vectorized in this fashion have 
a bright future. Because these architectures are iB general regular, the 
algorithms that cap be vectonzed for such machines tend to be regular,^ 

^i.e., simple, nonadapti^^, uniform-^d, low-order algorithms. It is clear 
thatjthere are nicely behaved problems for which these regular algofithms 
on ±>peciaUzed machines will require significantly less time than algorithms 
requiring fewer operations on serial machine. Il is also clear that no 
ipa^er how fest tt^ specialized machijpes ar^^ there are problems that 
are^fflcieflfl} diffiollt that more sophisticated algorithms are needed for 
more^^jSyal-purpose computers. * ^ 

In order to bring about these advances in architecture, it ts necessary t 
to involve practitioners of scientific computation in the design process. 
Luckily, modem design automationi tools should mUce it possible for in- 
terdisciplinaiy^ design teams to succeasfutly synthesize innovative special- 
purpose' systems, and automated fabrication facilities should make it po$^ 

' sible for such systems to hte built,, debugged, and used. ^ ' ^ * 

- These advances^ in macb ine a rchitecture should also have a dramatic ^ 

effect on the 'design an*i] analysis ofaTgoriCIims foracientific-computing. ' ^ 
TpaditionallV, such analysis ^as beenJ3ased on (asymptotic) estimates of _ * 
the number of arithmetic operationil^However, with these jew architec- 
tures it is quite likely that the ninnm^^time of an algorithm will be ipore 

ide^endent on the movement of data than on the number of arithmetic 

. Y)perations. Thus, we need to develop new analytic modelsof complexit}rof 
scientific algorithms so thatWch models give us usefj^l information about 
the relative performance of algorithms. 
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